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ABSTRACT 


The  goal  of  this  work  is  to  develop  methodologies  for  synchronizing  the  clocks  of 
neighboring  nodes  of  an  underwater  acoustic  network.  Clock  synchronization  requires  an 
estimate  of  the  signaling  time  delay  between  the  nodes,  which  is  particularly  challenging 
in  the  presence  of  multipath  propagation  through  the  acoustic  communication  channel. 
This  thesis  focuses  on  modeling  the  underwater  acoustic  communication  channel, 
accounting  for  the  multipath  arrivals,  and  creating  a  set  of  signal  processing  algorithms 
for  estimating  the  required  delay  times  that  enable  clock  synchronization  protocols  for  the 
underwater  acoustic  network.  The  proposed  method  involves  correlating  the  responses  of 
the  bidirectional  channels  to  exploit  the  underlying  reciprocity.  Performed  in  two  stages, 
a  sequence  of  probe  signals  is  first  transmitted  to  create  an  ensemble,  which  contains 
information  about  the  time-variability  of  the  acoustic  communication  channel  with 
multipath.  From  this  ensemble,  we  determine  its  dominant  time-invariant  characteristic 
and  use  it  as  a  reference  datum  for  the  time  delay  measurements.  The  second  stage 
consists  in  performing  time-delay  estimation  of  two  probe  signals  exchanged  between 
nodes.  The  two  stages  are  tested  using  simulated  signal  measurements,  and  actual  signal 
measurements  were  performed  in  a  fresh-water  lake  for  the  first  phase  only.  Both 
computer  simulations  and  experimental  results  show  the  effectiveness  of  the  proposed 
techniques. 
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EXECUTIVE  SUMMARY 


Clock  synchronization  among  network  nodes  becomes  necessary  when  data  from 
multiple  nodes  are  to  be  co-processed.  In  the  context  of  surveillance  networks, 
synchronization  of  the  nodes’  clocks  is  directly  related  to  the  accuracy  of  the  localization 
of  a  potential  target.  Clock  synchronization  is  achieved  by  evaluating  and  then  correcting 
the  skew  and  offset  of  the  internal  clock  of  a  single  node.  The  skew  of  a  clock  is  the 
oscillator  frequency  difference  between  two  synchronizing  clocks,  whereas  the  offset  is 
the  time  difference  at  a  specific  point  in  time  between  two  clocks. 

Clock  synchronization  is  especially  important  for  high-latency  communications 
such  as  those  in  underwater  acoustic  networks  in  which  data  from  separate  nodes  take 
variable  amounts  of  time  to  reach  the  main  processing  node.  The  challenge  in 
underwater  acoustic  clock  synchronization  lies  in  correctly  estimating  the  time  delay 
between  the  transmission  of  a  signal  from  a  node  and  the  reception  of  that  signal  by 
another  node.  For  the  acoustic  communications  problem,  the  propagation  medium 
introduces  environmental  factors  such  as  high  latency,  multipath,  scattering,  refraction, 
transmission  loss,  noise  and  non-reciprocity.  These  factors  have  significant  effects  on  the 
channel  impulse  response.  The  individual  multipath  arrival  amplitudes  and  the  arrival 
times  can  also  vary  greatly  depending  on  the  surface  conditions  and  the  directivity  of  the 
transducers.  For  some  environmental  conditions,  the  first  or  other  multipath  arrivals  may 
be  greatly  attenuated  or  absent  altogether. 

The  aim  of  this  thesis  is  to  examine  how  time-variability  of  the  impulse  response 
impacts  clock-synchronization  and  how  it  can  be  avoided  or  mitigated  in  a  spatially  static 
underwater  network.  The  Time  Synchronization  for  High  Latency  (TSHL)  protocol 
developed  by  Syed  and  Heidemann  in  [1]  is  used  as  the  basis  for  this  research.  The 
protocol  is  adapted  to  account  for  time-variability  of  the  channel  impulse  response. 

In  the  absence  of  a  fit-for-all  underwater  statistical  model,  channel  soundings  are 
a  good  alternative  to  the  lack  of  prediction  capability.  They  are  direct  measurements  of 
the  time-varying  impulse  response  of  a  propagation  medium  and  replace  general 


assumptions  about  the  prevailing  channel  between  two  nodes.  A  series  of  pulses  is  sent 
through  the  channel  at  a  set  interval.  The  receiver  analyzes  the  received  series  and 
compiles  its  own  local  channel  characteristics  from  it. 

Our  approach  relaxes  the  assumption  in  TSHL  of  a  stable  propagation  delay 
between  the  synchronizing  nodes.  It  also  relaxes  the  reciprocity  assumption  for  the 
channel.  Instead,  we  repeatedly  measure  the  response  of  the  channel  over  a  period  of 
time  to  statistically  define  the  communication  channel  to  better  estimate  the  propagation 
delay  and  its  skew.  This  approach  could  also  be  defined  as  adding  channel  sounding  to 
the  TSHL  protocol.  Our  model  keeps  the  two-stage  approach  developed  by  TSHL, 
correcting  for  the  skew  first  followed  by  the  offset. 

In  the  skew-correction  stage,  we  analyze  a  series  of  chirp  responses  to  identify  the 
least  time-varying  multipath  present  in  the  channel  between  the  two  nodes.  Based  on  the 
known  chirp  transmission  interval,  linear  regression  is  performed  on  the  measured 
interval  of  arrivals  of  that  identified  multipath  detected  within  each  chirp  response.  The 
slope  of  the  resulting  best-of-fit  straight  line  represents  the  skew  of  the  clock’s  receiver 
node  compared  to  the  originating  or  reference  node. 

The  second  stage  consists  of  an  exchange  of  time-stamp  messages  between  the 
two  nodes,  where  we  assume  a  time-varying,  one-way  propagation  delay  between  the  two 
nodes.  The  challenge  in  this  stage  consists  in  determining  the  difference  in  propagation 
delays.  This  is  accomplished  by  performing,  at  one  of  the  nodes,  a  time-delay  estimate 
using  the  chirp  responses  of  the  channel  sounding  of  both  one-way  propagation  paths. 

The  Bellhop  and  VirTEX  underwater  environment  models  are  used  to  simulate  a 
simple  communication  channel  in  Matlab.  Variability  is  introduced  in  the  model  by 
modifying  the  state  of  the  surface  and  rotating  one  of  the  nodes  between  the 
transmissions  of  the  chirps.  The  aim  of  the  simulations  was  to  prove  our  conceptual 
approach  prior  to  conducting  an  in-water  experiment  and  to  demonstrate  the  effect  of 
time-variability  on  our  clock-synchronization  process.  The  simulation  verified  the  theory 
of  our  approach  as  it  precisely  determined  the  difference  in  propagation  delays  and  clock 
offset. 
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The  first  stage  of  the  algorithm  was  then  tested  in  the  Del  Monte  Lake  situated  on 
the  NPS  campus.  The  general  concept  of  the  experiment  consisted  of  transmitting  and 
recording  a  series  of  chirps  between  two  acoustic  modems.  The  characteristics  (duration 
and  modulation)  of  the  chirp  were  modified  between  the  transmissions  of  the  transmitted 
chirp  series  to  analyze  the  impact  of  the  chirp  selection  on  the  protocol.  The  digital 
recordings  were  recovered  and  analyzed  using  the  same  algorithms  exercised  in  the 
simulation.  The  second  stage  of  the  algorithm  was  not  tested  since  it  would  require 
modifying  the  firmware  of  the  acoustic  modems  to  program  the  message  exchange 
protocol. 

Overall,  the  results  from  the  experiment  demonstrate  some  time  variability  in  the 
tested  channel.  The  first  arrival  is  always  present  and  strong.  However,  the  following 
arrivals  are  observed  to  be  gradually  less  stable  as  they  increasingly  interact  with  the 
channel  boundaries.  The  algorithm  is  able  to  identify  various  multipaths  based  on  the 
detected  arrivals  and  detennines  the  most  stable  one  based  on  the  correlation  coefficient 
of  a  model  fit  to  the  time-of-arrival  estimates. 

Additional  experiments  are  required  to  fully  prove  the  developed  algorithms  and 
select  the  optimal  parameters,  such  as  the  probe  signal  characteristics.  Once  completely 
tested,  the  addition  of  this  capability  to  an  underwater  acoustic  network  should  allow  it  to 
keep  the  clocks  of  the  nodes  synchronized.  This  synchronization  should  improve  the 
reliability  of  the  network  when  deployed  or  used  for  long  periods  of  time. 

[1]  A. A  Syed  and  J.  Heidemann,  “Time  Synchronization  for  High  Latency 

Acoustic  Networks,”  in  Proc.  Of  IEEE  Conference  on  Computer  Communications 
( Infocom ),  Barcelona,  Spain,  Apr.  2006. 
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I.  INTRODUCTION 


A.  REQUIREMENT  FOR  CLOCK  SYNCHRONIZATION 

Clock  synchronization  in  distributed  sensor  networks  has  been  a  research  topic  for 
many  years.  Synchronization  among  network  nodes  becomes  necessary  when  data  from 
multiple  nodes  are  to  be  co-processed.  In  the  context  of  target  detection  and 
triangulation,  the  synchronization  of  the  nodes’  clocks  is  directly  related  to  the  accuracy 
of  the  localization  of  a  potential  target. 

The  majority  of  recent  developments  in  the  literature  involve  wireless  radio 
networks  where  signal  propagation  is  almost  instantaneous  (speed  of  light).  In  most  of 
the  literature,  clock  synchronization  is  accomplished  via  the  exchange  of  time-stamped 
messages  between  two  nodes.  This  is  achieved  by  evaluating  and  then  correcting  the 
skew  and  offset  of  the  internal  clock  of  a  single  node.  The  accuracy  of  clock 
synchronization  is  directly  related  to  the  accuracy  of  the  measured  propagation  time 
delays  among  the  nodes. 

B.  CLOCK  SYNCHRONIZATION  IN  THE  UNDERWATER  ACOUSTIC 

CHANNEL 

Clock  synchronization  is  especially  important  for  high-latency  communications 
such  as  those  in  underwater  acoustic  networks  where  data  from  separate  nodes  take 
variable  amounts  of  time  to  reach  the  main  processing  node.  The  challenge  in 
underwater  acoustic  clock  synchronization  lies  in  correctly  estimating  the  time  delay 
between  the  transmission  of  a  signal  from  a  node  and  the  reception  of  that  signal  by 
another  node.  For  the  acoustic  communications  problem,  the  propagation  medium 
introduces  environmental  factors  such  as  high  latency,  multipath,  scattering,  refraction, 
transmission  loss,  noise,  and  non-reciprocity.  These  factors  have  significant  effects  on 
the  channel  response.  Prior  to  perfonning  clock  synchronization,  we  must  address  the 
complex  time  delay  estimation  issues  associated  with  the  response  of  the  channel.  The 
idealized  model  for  the  channel  impulse  response  is  a  sequence  of  impulses  representing 
the  multipath  arrivals  caused  by  the  channel.  In  theory,  we  can  estimate  the  time  delay 
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between  transmission  and  reception  by  identifying  a  reference  arrival  (e.g.  the  first  arrival 
or  another  reliable  arrival).  However,  the  individual  multipath  arrival  amplitudes  and  the 
arrival  times  can  vary  greatly  depending  on  the  surface  conditions  and  the  directivity  of 
the  transducers.  For  some  environmental  conditions,  the  first  or  other  multipath  arrivals 
may  be  greatly  attenuated  or  absent  altogether.  This  makes  the  identification  of  a 
reference  arrival  very  difficult. 

While  the  previously  stated  difficulties  associated  with  clock-synchronization 
protocols  in  an  underwater  network  are  important,  the  extreme  bandwidth  and  latency 
limitations  of  underwater  communications  constrain  the  implementation  of 
synchronization  protocols  that  involve  negotiation  or  arbitration.  Most  underwater 
acoustic  networks  can  only  achieve  very  low  data  communication  rates  and  care  must  be 
taken  to  limit  the  amount  of  information  exchanged  between  the  nodes  of  the  network. 
Associated  with  this  limitation  is  the  finite  power  supply  at  each  node.  To  maximize  the 
operational  utilization  of  a  node,  the  number  of  transmissions  and  their  associated 
message  length  must  be  limited.  In  the  context  of  clock  synchronization,  the 
communications  exchanged  between  the  synchronizing  nodes  must  be  minimized,  while 
the  length  of  time  between  the  need  to  repeat  the  clock  synchronization  algorithm  should 
be  maximized. 

The  aim  of  this  thesis  is  to  examine  how  time  variability  of  the  impulse  response 
impacts  clock-synchronization  and  how  it  can  be  mitigated  in  a  spatially  static 
underwater  network.  The  Time  Synchronization  for  High  Latency  (TSHL)  protocol 
developed  by  Syed  and  Heidemann  in  [1]  is  used  as  the  basis  for  this  research.  It  is  often 
referred  to  as  the  source  of  other  work  in  the  literature,  like  this  research. 

Originally,  TSHL  is  designated  for  spatially  fixed  underwater  networks.  In  this 
research  it  is  adapted  to  account  for  time  variability  of  the  channel  impulse  response  and 
tested  using  computer  simulations.  Actual  acoustic  signal  measurements  acquired  under 
controlled  conditions  in  a  fresh-water  lake  serve  to  verify  the  physical  assumptions  that 
are  the  basis  of  the  protocol. 
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This  thesis  is  organized  as  follows.  The  general  concepts  associated  with  clock- 
synchronization  and  the  variability  of  the  underwater  communication  channel  are 
introduced  in  Chapter  II  and  III.  Our  proposed  clock-synchronization  algorithms 
accounting  for  time -variability  in  the  channel  are  presented  in  Chapter  IV.  These 
algorithms  are  then  simulated  in  Chapter  V  and  some  of  them  tested  in  an  actual 
underwater  environment  in  Chapter  VI.  Finally,  the  recommendation  and  conclusions  are 
given  in  Chapter  VII  and  VIII. 
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II.  CLOCK  SYNCHRONIZATION 


A.  THE  COMPUTER  CLOCK 

A  computer  clock  is  an  ensemble  of  hardware  and  software  components  used  to 
provide  an  accurate,  stable  and  reliable  time-of-day  function  to  the  operating  system  and 
its  clients  [2].  From  a  hardware  point  of  view,  most  clocks  are  based  on  a  mechanism 
that  oscillates  at  a  precise  constant  frequency.  The  precision  of  the  oscillator  to  generate 
this  frequency  is  defined  as  its  tolerance,  and  since  it  is  usually  very  low,  it  is  expressed 
in  terms  of  parts-per-million  (ppm).  How  well  it  can  maintain  this  set  frequency  is 
defined  as  the  stability  of  the  clock.  The  oscillating  mechanism  can  take  the  fonn  of  a 
simple  mechanical  pendulum,  or  the  vibrations  of  electrons  in  an  atom.  These  two 
technologies  are  situated  at  opposite  ends  of  the  spectrum  in  terms  of  precision,  stability, 
complexity  and  cost.  In  the  vast  majority  of  cases,  the  middle-of-the-road  solution  is  to 
use  a  man-made  quartz  crystal  as  the  oscillating  mechanism. 

By  knowing  the  oscillating  frequency  of  the  mechanism,  a  unit  of  time  can  be 
linked  to  a  set  number  of  oscillation  periods.  In  the  case  of  a  computer  clock,  a  simple 
counter  will  send  an  interrupt  message  incrementing  the  clock  by  one  time  unit  following 
a  pre-determined  number  of  oscillations.  Sending  this  interrupt  resets  the  counter  to  0 
and  reinitiates  the  process. 

B.  THE  QUARTZ  CRYSTAL  OSCILLATOR 

1.  General  Description 

More  than  2  billion  quartz  crystal  oscillators  are  manufactured  annually  [2].  They 
are  used  in  consumer  electronics  such  as  wristwatches,  mobile  phones  and  computers, 
and  in  equipment  requiring  greater  precision  such  as  test  and  measurement  instruments. 
Quartz  crystal  is  abundant  in  nature  and  is  easy  to  grow  in  large  quantities,  at  low  cost, 
and  with  relatively  high  purity.  It  is  the  piezoelectric  property  of  the  quartz  crystal  which 
makes  it  a  material  of  choice  for  use  as  a  stable  resonator  in  oscillators.  Its  high  Q  value 
makes  it  an  excellent  material  for  use  as  an  oscillator  [2].  A  quartz  crystal  can  be 

manufactured  in  various  structures  and  cuts,  which  define  its  own  operational 
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characteristics.  One  such  characteristic  is  the  resonant  frequency  which  can  be  from  a 
few  tens  of  kilohertz  to  tens  of  megahertz. 

2.  Quartz  Crystal  Instabilities 

A  quartz  crystal  often  includes  impurities  which  lead  to  variations  in  the 
resonance  frequencies.  This  is  especially  true  for  less  expensive  oscillators  found  in  most 
electronic  devices,  where  the  frequency  can  vary  up  to  40  ppm  between  two  clocks  [3]. 
The  resonating  frequency  of  a  quartz  crystal  is  subject  to  many  types  of  instabilities 
which  can  affect  its  precision  and  accuracy  [4].  These  include  aging,  noise,  and 
frequency  changes  depending  on  temperature,  acceleration,  ionizing  radiation,  power 
supply  voltage,  etc.  By  knowing  how  these  instabilities  affect  a  clock,  it  is  possible  to 
correct  for  them  using  hardware  or  software  solution.  As  an  example,  the  real-time  clock 
used  on  Benthos  Teledyne  acoustic  modems  is  temperature  compensated  [5]. 

C.  THE  NEED  FOR  A  REFERENCE  TIME 

Every  single  device  in  a  network  will  have  at  least  one  local  clock,  the  time  of 
which  will  be  different  from  another  unless  we  provide  proper  synchronization. 
Depending  on  how  much  the  times  differ,  the  following  fundamental  operations 
performed  by  a  node  in  a  wireless  sensor  network  may  or  may  not  be  possible  [3]. 

There  are  a  number  of  issues  strictly  related  to  clock  management,  and  they  are 
listed  next. 

1.  Data  Fusion 

This  is  one  of  the  principal  benefits  to  be  derived  from  a  sensor  network,  where 
the  goal  is  to  observe  and  report  on  events.  To  coherently  integrate  and  process  event 
information  from  multiple  sensors  in  a  meaningful  way,  all  data  must  be  referenced  to  a 
standard  time. 

2.  Power  Management 

One  way  to  manage  and  save  energy  in  a  network  is  by  the  implementation  of  a 
duty  cycling  protocol.  Instead  of  leaving  a  node  operational  for  a  continuous  period  of 
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time,  it  can  be  programmed  to  sleep  and  wake-up  at  a  set  time  interval.  While  in  sleep 
mode,  the  node  saves  energy  by  keeping  only  essential  systems  operating  and  expending 
minimum  power.  To  permit  communication  between  two  nodes,  both  of  them  need  to  be 
awake  at  the  same  time.  Here  lies  the  importance  that  both  nodes  operate  on  a 
synchronized  clock,  such  that  they  can  simultaneously  switch  between  operating  modes. 

3.  Transmission  Scheduling 

Some  network  protocols  allow  a  node  access  to  a  network  by  implementing  a  set 
schedule  for  when  a  specific  node  can  transmit  (i.e.,  the  time-division  multiple  access 
protocol )  [3],  From  this  schedule,  only  one  node  can  transmit  at  a  time,  reducing  the  risk 
of  collisions  within  the  network.  Another  advantage  of  such  a  protocol  is  the  reduction  of 
the  message  overhead  that  is  usually  required  to  control  the  communication  flow  within 
the  channel,  therefore  increasing  the  available  bandwidth.  The  main  requirement  for  this 
kind  of  protocol  is  that  every  node  in  the  network  must  use  synchronized  clocks  to  start 
and  stop  their  transmissions  on  a  common  distributed  schedule. 

4.  Master-slave  Relationship 

The  above  operations  have  a  mutual  requirement:  the  internal  reference  clock  of 
every  node  in  the  network  must  be  synchronized.  For  synchronization  to  exist,  the 
concept  of  a  master-slave  relationship  must  be  introduced  in  the  network.  In  the  case  of  a 
sensor  network,  a  node  is  designated  as  master  and  its  internal  clock  is  used  as  a 
reference  by  all  the  other  nodes,  designated  as  slaves.  To  avoid  having  to  repeatedly  re¬ 
synchronize  the  clocks  of  the  network,  the  master  node  is  often  equipped  with  a  more 
stable  and  precise  clock  than  the  other  nodes  or  will  have  the  capability  to  synchronize 
with  an  external  reference  time  such  as  the  one  provided  by  the  GPS.  Having  an  accurate 
master  clock  is  desirable  in  a  sensor  network  but  is  not  as  critical  as  having  all  the  clocks 
of  the  network  synchronized  with  the  master. 
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D. 


CLOCK  MODEL 


Each  device  or  node  in  a  sensor  network  has  its  own  clock.  The  time  measured  by 
clock  A  usually  presents  random  errors,  but  on  average  it  can  be  represented  by  a  linear 
function  of  the  form 

CA(t)  =  fAt  +  bA  (1) 

where  the  temi  fA  represents  the  skew  and  bA  the  offset  of  clock  CA  compared  to  its 
reference  or  master  clock  CB  .  Throughout  this  paper  B  denotes  the  reference  node  and  A 
the  synchronizing  node.  In  an  ideal  situation  with  no  errors  or  drifts,  we  have  C(t)  =  t . 
The  skew  of  a  clock  fA  is  the  oscillator  frequency  difference  between  CA  and  CB  .  The 
offset  bA  is  the  time  difference  at  a  specific  point  in  time  between  the  two  clocks.  Figure 
1  is  a  graphical  representation  of  a  perfect  reference  clock  combined  with  both  a  slower 
and  faster  slave  with  respective  offset  clocks. 


Figure  1.  Clock  model  of  a  fast  (C2),  reference  (Ci),  slow  clock  (C3)  and  their 

associated  equations. 

E.  SOURCES  OF  ERROR  IN  CLOCK  SYNCHRONIZATION 

In  a  perfect,  ideal  environment  with  no  delayed  responses,  clock  synchronization 
is  simple.  Node  B  sends  its  current  time  to  node  A,  which  adjusts  its  own  clock 
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accordingly.  In  an  actual  network,  many  steps  are  required  to  accomplish  this  exchange 
of  information.  Associated  with  each  of  these  steps  is  a  systemic  time  delay  which  can 
be  either  detenninistic  or  stochastic. 

The  following  is  a  list  of  the  various  delays  generated  by  the  transmission  of  a 
timing  message  between  two  nodes  [1]  and  illustrated  using  the  open  systems 
interconnection  (OSI)  model  in  Figure  2.  The  OSI  model  structures  a  communication 
system  in  terms  of  functional  layers.  Each  layer  can  only  communicate  with  the  layer 
above  or  under  it  within  a  communication  device.  In  order  for  a  communication  system 
to  exist,  two  devices  or  nodes  must  share  the  same  layered  structure. 

1.  Send  Time 

Send  time  is  the  time  delay  from  generating  the  timing  message  at  the 
synchronizing  node  by  its  software  as  transfers  through  the  application,  transport  and 
network  layers  on  the  sender  side  illustrated  in  Figure  2.  This  delay  results  from  the 
execution  of  software  operations  and  is  a  fixed  time  in  most  devices. 

2.  Access  Time 

This  is  the  contention  time  to  access  the  channel.  It  is  highly  variable  and 
depends  on  the  MAC  protocol  used  and  how  busy  the  channel  is  when  the  node  is 
attempting  to  access  it. 

3.  Transmission  and  Reception  Time 

This  is  the  time  to  send  and  receive  the  entire  message  by  the  transducer. 
This  is  a  detenninistic  delay,  which  depends  on  the  size  of  the  message  sent  and  the 
bandwidth  of  the  link.  This  is  illustrated  by  the  four  green  squares  in  Figure  2,  each 
representing  a  packet. 
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4.  Propagation  Time 


The  propagation  is  the  time  required  for  the  timing  message  to  travel 
between  the  two  nodes.  The  high-latency  nature  of  water  as  a  communication  medium 
makes  detennining  this  delay  very  difficult.  This  problem  is  further  explained  in  Chapter 
III. 


5.  Encoding  and  Decoding  Time 

This  is  the  time  required  for  the  digital  signal  processor  (DSP)  in  the 
acoustic  modem  to  encode/decode  the  message  to/from  the  received  acoustic  signal.  This 
is  a  short  and  detenninistic  delay. 

6.  Receive  Time 

Similar  to  the  send  time,  this  is  the  delay  for  the  software  to  process  the 
received  timing  message  through  the  upper  layers  of  the  OSI  model  as  illustrated  in 
Figure  2. 
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Figure  2.  Sources  of  error  in  estimating  message  latency  explained  using  the  OSI 

model.  After  [1]. 


To  measure  and  eliminate  some  of  the  uncertainties  associated  with  these  delays, 
multiple  messages  need  to  be  exchanged  between  the  two  nodes,  adding  yet  more  delays 
and  complexity  to  the  overall  process.  In  some  sense,  clock  synchronization  in  a  sensor 
network  can  be  regarded  as  the  process  of  accounting  for  the  effects  of  unknown  delays 
in  the  overall  communication  process  [3]. 

F.  CLOCK  SYNCHRONIZATION  BASICS 

From  Equation  1,  we  see  that  two  parameters  (skew  and  offset)  need  to  be 
corrected  to  synchronize  a  local  clock  to  a  reference.  In  general,  the  process  to  execute 
this  correction  is  achieved  by  transferring  a  set  of  timing  messages  between  the  two 
nodes  as  shown  in  Figure  3  and  referred  to  as  the  two-way  message  exchange  (or  sender- 
receiver  synchronization)  [6]. 
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Synchronizing 
Node  A 


Reference 
Node  B 


Figure  3.  Two-way  message  exchange  of  timing  messages  between  the  synchronizing 

(A)  and  reference  (B)  nodes.  After  [1], 


Most  authors  use  Figure  3  to  describe  the  basics  of  clock  synchronization 
protocols.  The  literature  usually  assumes  a  reciprocal  and  symmetric  communication 
channel  where  the  one-way  propagation  delay  can  be  determined  by  averaging  the  total 
propagation  time  of  the  two  message  exchanges.  Why  Figure  3  can  be  amended  to 
illustrate  a  time -varying  propagation  delay,  such  that  the  signal  propagation  delay 
between  node  A  and  B  can  be  different  from  the  one  from  B  to  A,  is  discussed  in  the  next 
chapter.  A  solution  to  address  this  non-reciprocity  is  proposed  in  Chapter  IV. 

From  Figure  3,  the  total  propagation  time  can  be  mathematically  modeled  as 

(T4-Tl)  —  (T3-T2)  =  2p  =  P.  (2) 

Assuming  a  constant  propagation  delay  and  zero  or  negligible  skew,  we  can 
detennine  the  offset  of  CA  from  a  single  exchange  of  time-stamped  messages  such  as 
illustrated  in  Figure  3.  Once  the  propagation  delays  have  been  determined,  the  offset  bA 
can  be  calculated  from  one  of  the  following: 


T2-T\-  p  =  bA 

(3) 

T4-T3-  p  -  —bA . 

(4) 
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The  offset  bA  from  Equation  3  changes  sign  in  Equation  4  since  the  offset  sign  is 
relative  to  which  node  it  is  compared  with. 

By  repeating  the  exchange  of  time-stamped  messages,  the  clock  synchronization 
protocol  can  also  determine  the  clock  skew  using  this  fonnulation.  The  Time 
Synchronization  for  Time-Varying  Propagation  Underwater  Network  (TSVP)  [7]  keeps 
track  of  the  time-stamps  and  calculates  the  clock  skew  using  linear  regression. 

Another  well-known  timing  signaling  approach  is  the  one-way  message 
dissemination  method  [6],  The  master  node  broadcasts  its  timing  infonnation  to  all  the 
nodes  in  its  vicinity  at  a  set  time  interval.  A  slave  node  collects  and  saves  all  the  timing 
information  in  a  matrix  to  calculate,  using  a  least-squares  estimator,  the  skew  and  offset 
of  the  local  clock  compared  to  the  reference  clock.  This  method  assumes  a  fixed 
propagation  delay,  but  the  effect  of  this  assumption  is  mitigated  by  the  fact  that  the  least- 
squares  estimator  in  fact  averages  these  delays.  This  methodology  is  used  by  the 
flooding  time  synchronization  protocol  (FTSP)  [8],  Although  this  approach  requires 
multiple  broadcast  messages  from  the  master  node,  it  can  synchronize  a  multitude  of 
slave  nodes  simultaneously  and  requires  minimal  channel  access  negotiation  between  the 
nodes. 

How  these  two  signaling  approaches  are  combined  and  used  in  two  different 
stages  to  produce  a  clock  synchronization  protocol  is  described  in  Chapter  IV. 

G.  LIMITATION  TO  CLOCK  SYNCHRONIZATION 

It  is  possible  to  keep  all  local  clocks  of  a  network  very  precisely  synchronized  by 
simply  updating  them  on  a  periodic  basis.  However,  this  scenario  is  unrealistic  and 
would  compromise  the  efficient  use  of  the  node’s  two  limited  resources  [6],  energy  and 
bandwidth,  as  explained  in  the  following. 

1.  Energy 

In  general,  devices  found  in  sensor  networks,  even  underwater  ones,  are  limited  in 
size  and,  therefore,  restricted  in  the  amount  of  energy  they  can  hold.  In  order  to 
maximize  its  battery  and  operational  life,  clock  synchronization  must  be  achieved  while 
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preserving  energy  for  normal  operations.  Since  transmitting  is  an  energy-consuming 
operation,  a  clock-synchronization  protocol  must  limit  the  number  and  the  length  of  the 
messages  exchanged  within  the  network. 

2.  Bandwidth 

Underwater  networks  have  limited  bandwidth  and  are  restricted  to  a  relatively 
slow  transmission  data  rate.  Clock  synchronization  is  a  network  support  feature  and  not 
the  main  role  of  a  network  node.  Limiting  the  bandwidth  used  for  clock  synchronization 
helps  to  maximize  the  utilization  of  the  node  for  its  intended  main  purpose. 

A  balance  based  on  the  trade-offs  between  the  desired  clock  synchronization 
precision  and  the  limited  network  resources  must  ultimately  be  achieved  to  maximize  the 
efficiency  of  the  network. 
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III.  VARIABILITY  OF  THE  UNDERWATER  COMMUNICATION 

CHANNEL 


The  underwater  acoustic  channel  is  a  challenging  communication  medium.  It  can 
be  thought  of  as  a  long  slim  pipe,  with  the  poor  physical  link  quality  of  a  mobile 
terrestrial  radio  channel,  and  the  high  latency  of  a  satellite  channel  [9],  Technologies 
developed  for  wireless  radio  communications  are  often  used  as  the  basis  for  underwater 
acoustic  communications.  The  ocean  shares  with  the  over-the-air  medium  many  of  the 
same  impairments  that  make  wireless  communications  difficult  such  as  multipath, 
scattering,  noise,  attenuation,  and  fading.  The  difficulties  introduced  by  these 
characteristics  are  magnified  in  the  ocean  due  to  the  high  latency  of  the  medium. 
Acoustic  signals  propagate  five-orders-of-magnitude  slower  in  water  compared  to  over- 
the-air  RF  signals,  with  speeds  of  1.5xl03  m/s  and  3xl08  m/s,  respectively.  In  addition, 
greater  absorption  at  higher  acoustic  frequencies  limits  the  available  bandwidth,  leaving 
only  lower  frequencies,  which  are  susceptible  to  channel  fading  due  to  multipath. 

Underwater  acoustic  research  usually  takes  into  account  the  above-mentioned 
intricacies  of  the  medium.  However,  most  practical  solutions  make  the  assumption  that 
the  environment  is  stable  over  relatively  short  periods  of  time.  Nevertheless,  signal 
fluctuations  can  occur  due  to  transceiver  motion  or  inherent  changes  within  the 
propagation  medium.  These  changes  are  due  to  numerous  phenomena  which,  depending 
on  their  time  scale,  might  or  might  not  affect  the  instantaneous  quality  of  the 
communication  channel.  They  will  also  have  a  different  impact  depending  on  whether 
the  channel  is  in  deep  or  shallow  water.  The  shallow-water  channel  is  recognized  as 
being  the  most  difficult  due  to  its  time-varying  impulse  response.  This  greatly  affects  the 
signal  propagation  delay  and  impacts  the  clock  synchronization  effort. 

A.  TRANSDUCER  MOTION 

In  the  underwater  environment  there  is  always  some  motion  present  in  the  system 
[9].  Motion  can  affect  the  impulse  response’s  frequency  spectrum  (via  the  Doppler 
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effect),  strength  of  the  arrival  paths,  or  in  the  worst  case,  the  complete  disappearance  of  it 
due  to  one  of  the  transceivers  entering  a  shadow  zone  where  no  communication  path 
exists  between  two  communicating  nodes. 

The  Doppler  effect  will  cause  frequency  shifting  as  well  as  additional  frequency 
spreading.  The  amount  of  spectral  shift  is  proportional  to  the  ratio  of  the  relative 
transmitter-receiver  velocity  to  the  speed  of  sound  [9].  Since  the  speed  of  sound  in  the 
water  is  relatively  slow  (1500  m/s),  a  transducer’s  motion  from  waves,  currents,  tides, 
and  platform  motion  will  be  enough  to  impair  the  communication  channel.  Furthennore, 
unlike  wireless  communications  systems  in  the  air,  it  has  been  shown  that  assumption 
that  the  underwater  acoustic  channels  are  subject  to  wide-sense  stationary  uncorrelated 
scattering  (WSSUS)  does  not  always  apply  [10]. 

The  rotation  of  a  transducer  is  another  dimension  of  motion  that  can  affect  the 
impulse  response  of  the  channel.  As  shown  in  Figure  4,  in  a  multipath  environment,  the 
beam  pattern  of  a  transducer  influences  the  directivity  index  for  both  launch  angle  and 
receive  angle.  The  strongest  multipath  arrival  may  be  at  the  start  of  the  impulse  response, 
at  the  end,  or  somewhere  in  between  [10].  This  effect  is  demonstrated  in  Figure  4,  which 
illustrates  a  simple  multipath  model  of  an  environment  and  its  associated  impulse 
response  for  two  different  transducer  orientations.  Rotating  the  transducer  upward  by  30 
degrees  changes  the  relative  magnitudes  of  the  arrivals  in  the  channel’s  impulse  response, 
resulting  in  the  strongest  arrival  changing  from  the  2nd  to  the  1st  multipath  arrival. 
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Beampattern  of  transducer  at  0  degree  Beampattern  of  transducer  at  30  degrees 
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Figure  4.  Variation  of  the  impulse  response  caused  by  rotating  the  transducer  (and 
beampattern).  From  top  left  to  bottom  right:  simulated  sound-speed  profile 
(isothennal),  simulated  multipath  of  the  environment,  beampattern  when 
transducer  is  at  0  degrees,  beampattern  when  transducer  is  at  30  degrees,  impulse 
response  at  0  degrees  and  impulse  response  when  at  30  degrees.  Note  that  at  zero 
degrees,  the  second  arrival  has  the  strongest  magnitude,  and  at  30  degrees,  the 
strongest  magnitude  is  at  the  first  arrival. 


B.  TEMPORAL  VARIABILITY 

In  addition  to  transducer  motion,  the  impulse  response  will  vary  due  to  motion 
within  the  environment.  There  are  various  phenomena  responsible  for  this  motion,  which 
fall  in  the  two  following  categories  [9]  [11]. 
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1.  Long-Term 

Changes  that  occur  over  a  long  time  scale  usually  develop  at  a  slow  rate  and  do 
not  affect  the  signal  during  its  propagation  through  the  medium.  Their  gradual 
occurrence  means  that  the  system  can  adapt  to  the  changes  in  the  environment.  Since 
some  of  them  can  be  cyclic  in  nature,  the  system  can  be  designed  to  predict  and  adapt  to 
these  fluctuations. 


a.  Swell 

The  effect  of  swell  is  dominantly  indirect,  as  it  causes  the  motion  of  the 
transducers  and  affects  the  impulse  response  in  many  ways  as  explained  in  Section  A  of 
this  chapter.  It  can  also  directly  affect  the  signal  reflected  at  the  surface,  inducing  a 
Doppler  effect  on  the  signal  frequencies  based  on  the  direction,  speed  and  amplitude  of 
the  swell. 


b.  Current 

Spatially,  currents  change  the  sound-speed  profde  of  the  water  column. 
Temporally,  the  change  in  speed  of  the  propagation  medium  intrinsically  affects  the 
speed  of  propagation  of  the  signal.  This  change  of  speed  induces  a  Doppler  effect  on  the 
signal  frequencies.  The  current  creates  “scintillations”  (amplitude  instabilities),  due  to 
the  emergence  of  local  turbulence  in  the  current  stream. 

c.  Internal  Waves 

These  waves  have  similar  effect  as  the  current.  Their  inherent  movement 
induces  a  Doppler  effect  on  the  signal  frequencies. 

d.  Tides 

Their  period  of  half  a  day  makes  the  change  of  tides  to  the  environment 
relatively  slow.  They  only  affect  communication  channels  located  in  or  near  a  shallow- 
water  area,  where  they  modify  the  geometry  of  the  multipath  propagation  within  the 
medium. 
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e.  Daily  and  Seasonal  Temperature  Variations 

Temperature  variations  such  as  the  development  of  a  warm  surface  layer 
affect  the  sound-speed  profile  of  the  water  column.  Again,  their  rate  of  change  is  very 
slow  and  does  not  affect  the  instantaneous  characteristics  of  the  communication  channel. 

2.  Short-Term 

Short-term  variability  is  expressed  in  many  different  ways,  which  is  the  reason 
there  is  no  consensus  on  a  statistical  characterization  of  the  underwater  acoustic  channel 
[9].  Wind  energy  is  the  dominant  cause  of  short-term  instability.  Its  highly  dynamic 
nature  leaves  little  time  for  the  system  to  adapt  and  can  often  result  in  the  outage  of  the 
communication  channel. 

a.  Surface  Wave 

The  ever-moving  shape  of  the  ocean  surface  consistently  changes  the 
geometry  of  the  surface  reflection,  affecting  the  delays  of  the  impulse  response.  From 
Snell’s  law,  each  wave  crest  forms  an  acoustic  mirror  with  its  own  characteristic  shape 
and  focusing  properties  [12].  As  shown  in  Figure  5,  this  change  in  geometry  can  also 
create  caustics  resulting  in  signal  amplitudes  much  greater  than  those  of  the  direct  arrival 
path.  The  resulting  impulse  response  can  vary  in  time  in  such  a  manner  that  there  is  no 
clearly  defined  strongest  arrival  [9]. 


Figure  5.  Qualitative  description  of  acoustic  focusing  by  surface  wave  (from  [12]). 
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The  movement  of  the  wave  will  also,  as  in  the  case  of  swells,  imprint  a 
Doppler  effect  on  the  signal  frequencies. 

b.  Bubbles 

Air  bubbles  created  and  suspended  at  the  air-water  interface  result  in  a 
significant  attenuation  of  surface  scattered  signals,  especially  at  higher  signal  frequencies. 
The  “roughness”  of  the  sea-surface  breaking  waves  is  directly  related  to  the  amount  of 
bubbles  trapped  in  the  water  [13]. 

c.  Others 

There  is  a  long  list  of  diverse  occurrences  that  can  cause  time  variability 
and  affect  the  impulse  response  in  similar  ways  as  explained  previously.  Examples  of 
such  occurrence  are:  wakes  of  passing  ships,  fish  shoals,  bubble  screens  due  to  rain 
showers,  a  fresh- water  front  due  to  river  discharge,  etc  [10].  These  examples  are  often 
site-specific  and  neglected  when  designing  an  underwater  communication  system. 

C.  MEASUREMENT  NOISE 

Measurement  noise  is  a  problem  that  exists  in  and  affects  practically  every 
communication  channel.  However,  interference  due  to  noise  can  induce  errors  in  the 
reception  of  a  transmitted  signal  which  potentially  leads  to  errors  in  measuring 
propagation  time-delays  used  in  clock-synchronization  protocols.  Indeed,  successful 
communications  requires  a  link  margin  related  to  the  signal-to-noise  ratio  at  the  receiving 
node.  Like  any  sonar  system,  the  variability  of  clock  synchronization  can  be  analyzed 
with  the  sonar  equation. 

D.  CHANNEL  MODELING  AND  SOUNDING 

Effective  single  model  representations  of  the  salient  propagation  characteristics  of 
the  underwater  environment  have  been  elusive  [9].  Ocean  environment  modeling  is  still 
an  extremely  active  field  of  research  where  many  aspects  of  its  physics  are  still  not  fully 
understood.  Nevertheless,  the  recent  advances  in  computer  processing  power  have  made 
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it  possible  to  better  simulate  the  communication  channel.  These  simulations  are  based  on 
various  models,  which  are  not  always  valid  for  a  real  environment. 

In  the  absence  of  a  fit-for-all  underwater  statistical  model,  channel  soundings  are 
a  good  alternative  to  the  lack  of  prediction  capability.  They  are  measurements  of  the 
time-varying  impulse  response  of  a  propagation  medium  [10]  and  replace  general 
assumptions  about  the  channel.  A  series  of  pulses  is  sent  through  the  channel  at  a  set 
interval.  The  receiver  analyzes  the  received  series  and  estimates  the  local  channel 
characteristics  from  it.  This  procedure  is  an  additional  support  cost  to  the  system,  but  the 
result  may  greatly  improve  its  overall  efficiency. 
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IV.  PROPOSED  CLOCK  SYNCHRONIZATION  ALGORITHMS 


The  following  algorithms  account  for  time  variability  in  the  underwater 
communication  channel  when  synchronizing  the  clock  of  slave  node  A  to  reference 
master  node  B.  As  stated  in  the  introduction,  this  approach  is  based  on  the  TSHL 
protocol,  which  is  explained  in  the  following  section.  Our  approach  relaxes  the 
assumption  in  TSHL  of  a  stable  propagation  delay  between  the  synchronizing  nodes. 
Instead,  we  repeatedly  measure  the  response  of  the  channel  over  a  period  of  time  and 
statistically  describe  the  communication  channel  to  better  detennine  the  propagation 
delay  and  its  skew. 

Our  model  keeps  the  two-stage  approach  developed  by  TSHL,  but  the  stages  are 
implemented  differently  to  accommodate  for  their  respective  timing  message  signaling 
approaches  (one-way  message  dissemination  and  two-way  message  exchange).  The 
following  sections  describe  the  original  TSHL  protocol  in  more  detail  followed  by  an 
explanation  of  how  our  approach  extends  each  of  the  protocol’s  stages  to  account  for  a 
communication  channel  characterized  by  a  time-varying  response. 

A.  TSHL  ACOUSTIC  NETWORK  PROTOCOL 

The  TSHL  clock  synchronization  protocol  [1]  is  one  of  the  first  developed  and 
published  and  was  specifically  designed  for  underwater  or  high-latency  sensor  networks. 
Its  main  feature  is  to  split  the  clock  synchronization  process  into  two  distinct  stages  based 
on  the  importance  of  correcting  the  skew  of  a  clock  before  its  offset.  By  first  correcting 
the  skew  of  the  clock,  the  offset  is  converted  from  a  variable  function  of  time  to  a 
constant.  This  is  necessary  because  otherwise  the  high  latency  of  the  channel  makes  it 
impossible  to  precisely  correct  the  local  clock  since  its  offset  will  already  have  changed 
by  the  time  the  correction  is  applied. 

In  the  first  stage,  the  skew  of  the  clock  is  determined  using  a  one-way 
dissemination  method.  A  series  of  chirps  (x1,x2,xi,...,xn  )  is  sent  from  reference  node  B 
at  a  known  and  constant  reference  pulse  interval  (RPI)  to  the  synchronizing  node  A. 
Linear  regression  based  on  a  least-squares  fit  [14]  is  performed  on  the  measured 
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synchronization  pulse  interval  (SPI),  which  is  the  time  interval  between  the  arrival  time 
of  the  pulses  ( yt,y2,yi,...,yn),  relative  to  the  known  RPI.  The  resulting  slope  fA  of  this 

line  represents  the  clock  skew  of  node  A.  This  one-way  exchange  of  pulses  is  illustrated 
in  Figure  6.  The  first  stage  is  independent  of  the  propagation  delay  which  is  assumed  to 
be  constant  for  the  short  time  duration  when  the  series  of  chirps  is  transmitted. 


Figure  6.  One-way  dissemination  method  used  in  the  first  phase  of  the  TSHL  protocol 
to  correct  the  skew  of  the  local  clock  A  compared  to  the  reference  clock  B. 

When  stage  one  is  completed,  the  skew  fA  between  the  node’s  clocks  has  been 
determined,  and  bA  is  the  only  remaining  unknown  value  in  Equation  1.  The  known 
skew  can  be  corrected  such  that  clock  A’s  time-line,  illustrated  in  Figure  6,  is  now 
parallel  to  the  reference  clock  line  (horizontal).  At  this  point,  the  offset  bA  changes  from 
being  a  time-dependent  variable  to  a  constant  since  the  distance  between  the  two  lines 
remains  the  same. 

This  constant  offset  is  corrected  using  a  two-way  exchange  methodology.  As 
shown  in  Figure  3  and  assuming  a  constant  bi-directional  propagation  delay,  a  set  of 
messages  is  exchanged  from  node  A  to  B  and  then  from  B  to  A.  The  respective  node 
records  and  generates  time-stamps,  using  its  own  clock,  when  the  messages  are 
transmitted  ( 7,  and  T-, )  and  received  (T2  and  7’4  ).  This  timing  infonnation  is  compiled 


25 


by  the  synchronizing  node,  where  the  total  propagation  time  of  the  two  messages  will  be 
P  from  Equation  2.  Under  the  assumption  that  the  channel  is  reciprocal,  the  one-way 
propagation  time  delay  pk ,  where  k  is  the  message  number,  is  obtained  from  the  average 
of  the  two  messages 

=  (5) 

Knowing  the  one-way  propagation  delay  from  Equation  5,  where  pk-  p  ,  we  can 
use  Equations  3  or  4  to  determine  the  clock  offset  bA  .  Both  equations  should  yield  the 
same  value  for  the  offset  bA  . 

Eventually,  various  external  instabilities  (listed  in  Chapter  II)  will  affect  the 
oscillating  frequency  of  one  or  both  nodes,  changing  the  skew  difference  between  their 
local  system  clocks  thus  requiring  a  repeat  of  the  clock-synchronization  protocol. 

B.  STAGE  1:  CHANNEL  SOUNDING  AND  SKEW  CORRECTION 

Our  approach  to  the  skew  correction  looks  deeper  into  the  physical  attributes  of 
the  received  chirps.  When  subject  to  a  multipath  environment,  an  acoustic  modem  will 
“pick”  an  arrival  or  path  in  accordance  to  a  prc-dcfincd  metric.  In  the  case  of  the 
Teledyne  Benthos  acoustic  modems,  the  most  energetic  path  or  largest  peak  is  selected  as 
the  synchronizing  reference  datum.  As  explained  in  Chapter  III,  this  reference  arrival  can 
temporally  change  within  the  impulse  response  of  a  time-varying  communication 
channel.  Such  an  environment  will  influence  the  arrivals  within  the  impulse  response 
differently,  resulting  in  some  arrivals  being  more  stable  over  a  period  of  time  than  others. 
The  following  model  essentially  changes  the  prc-defined  synchronization  metric  used  by 
the  acoustic  modem.  The  new  metric  is  based  on  the  structure  of  the  impulse  response 
estimated  by  sending  a  series  of  chirps  through  the  channel  (channel  sounding). 
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Figure  7.  Skew  determination  system  diagram  illustrating  the  transmission  of  the  ith 

chirp  of  n  transmitted.  The  response  of  every  chirp  is  collected  and  recorded  by 
the  peak  detector  where  the  analysis  of  the  ensemble  begins,  resulting  in  the 
determination  of  the  skew  of  the  local  clock. 

The  block-diagram  of  the  proposed  system  is  illustrated  in  Figure  7. 

The  starting  point  of  the  process  is  the  generation  by  the  reference  node  B  of  a 
digital  impulse  train,  ultimately  based  on  its  clock  frequency.  The  length  of  this  impulse 
train  must  be  sufficient  to  provide  enough  chirps  to  the  synchronization  node  A  for  the 
linear  regression  algorithm.  Results  from  [1]  indicate  that  the  empirical  rule  of  thumb  for 
the  optimal  number  of  impulses  required  to  calculate  the  skew  is  25.  The  selection  of  the 
RPI  must  be  based  on  the  sampling  frequency  of  the  system  and  the  intended 
measurement  precision  of  the  skew.  The  length  of  the  whole  impulse  train  must  be  such 
that  there  are  enough  samples  to  cover  the  scale  of  the  clock  skew. 

The  characteristics  of  the  pulse  affect  how  accurately  each  multipath  arrival  can 
be  distinguished  from  one  another  within  the  impulse  response.  Since  the  function  of  the 
algorithm  is  essentially  the  same  as  channel  sounding,  reference  [10]  identifies  linear 
frequency-modulated  (LFM)  chirp  trains  and  pseudorandom  (PRN)  binary  sequences  as 
the  most  suitable  types  of  probe  signals.  A  key  attribute  of  LFM  and  PRN  wavefonns  is 
their  use  of  a  wide  spectral  bandwidth.  Another  factor  contributing  to  the  resolution  of 
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the  arrivals  is  the  length  of  each  chirp  by  virtue  of  the  large  time  and  bandwidth  product 
of  the  probe  signal.  A  longer  chirp  helps  identifying  the  pulse  within  a  noisy 
environment  but  depending  on  the  spread  of  the  impulse  response,  a  long  pulse  can  also 
introduce  some  superposition  in  the  sense  that  successive  pulses  in  time  cannot  be 
distinguished  from  each  other.  The  multipath  structure  and  distance  between  the  two 
synchronizing  nodes  determines  the  delay  spread  of  the  impulse  response  and  is, 
therefore,  linked  to  establishing  the  length  of  the  pulse  signal. 

As  the  series  of  chirps  is  created,  it  is  mixed  to  the  carrier  frequency  and 
transmitted  as  the  ith  chirp  ( x^t) )  of  n  total  chirps  by  the  reference  node’s  transducer. 
The  environment  affects  the  transmitted  signal  xi(t) ,  as  explained  in  Chapter  III.  The 
received  signal  is  then  denoted  as  yt(t)  in  Figure  7. 

The  signal  yt(t)  received  by  the  synchronizing  node’s  transducer  is  filtered, 
mixed  back  down  to  baseband,  and  digitized  based  on  the  sampling  frequency  regulated 
by  the  node’s  local  clock. 

The  time  delay  between  a  signal  recorded  at  two  spatially  separated  sensors  can 
be  recovered  using  a  correlator,  even  in  the  presence  of  noise  [15],  In  this  system,  the 
known  reference  pulse  xt(t)  is  correlated  with  the  signal  received  by  the  synchronizing 

node  yt(t)  to  produce  Rxy  (r) .  The  distance  between  the  peaks  in  Rxx  (r)  represent  the 

time  delays  between  the  multiple  copies  of  the  transmitted  pulse.  These  originate  from 
the  multipath  characteristic  of  the  channel.  (t)  can  also  be  used  to  estimate  the 

channel  impulse  response  at  time  t . 

The  peak  detector  then  applies  a  threshold  to  detect  the  location  of  the  peaks 
present  in  the  signal  R  (t).  These  peaks  represent  the  signal  propagation  time 

difference  between  the  different  paths  present  in  the  channel  at  time  t . 

Each  arrival  time  is  noted  in  Figure  7  as  J  which  refers  to  the  /h  received 

multipath  arrival  time  from  the  ith  out  of  n  transmitted  pulses.  The  arrival  times  are  then 
stored  in  a  matrix  A  as  illustrated  in  Figure  8.  It  is  important  to  note  that  the  number  of 
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detected  arrival  can  vary  among  the  received  chirp  responses.  Therefore,  m  represents 
the  maximum  number  of  arrivals  detected  in  a  single  received  pulse.  The  final  size  of 
matrix  A  is  then  n  x  m. 


j  =  multipath  arrival  times 
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Figure  8.  n  x  m  matrix  A  of  arrival  times  generated  by  the  peak  detector. 

The  peak  detector  will  not  necessarily  detect  the  same  number  of  arrivals  for  each 
received  pulse.  As  explained  in  Chapter  III,  the  time  variability  of  the  channel  impulse 
response  can  affect  the  energy  level  received  from  a  particular  multipath  and  prevent  it 
from  being  detected  by  the  peak  detector.  Therefore,  the  detected  arrival  times  present  in 
a  single  column  of  Figure  8  do  not  necessarily  represent  the  same  multipath  arrival  for 
each  one  of  the  received  responses.  The  arrivals  of  each  multipath  must,  therefore,  be 
properly  grouped  before  the  skew  can  be  found.  To  do  so,  a  sorting  algorithm  aligns  each 
multipath  in  a  single  column  of  a  matrix,  keeping  the  received  detected  arrival  for  a 
particular  pulse  on  the  same  row. 

To  detennine  the  main  multipath  arrivals,  a  probability  density  function  (PDF)  is 
estimated  using  all  of  the  detected  arrival  in  matrix  A  of  Figure  8.  The  time  indexes  of 
the  peaks  of  this  PDF  are  extracted  and  used  to  represent  the  main  multipath,  identified 
individually  by  the  index  q  in  Figure  9,  of  the  channel  where  q  =  1,2,. ..,<2 .  The 
maximum  number  of  multipaths  identified  by  the  algorithm  can  be  pre-set.  In  Figure  9, 
the  s  largest  peaks  of  the  PDF  will  be  identified  as  the  main  multipath,  where  s  =  Q  . 

Each  main  multipath  is  represented  by  a  column  in  a  new  matrix  B,  an  example  of 
which  is  illustrated  in  Figure  9.  The  time  indexes  are  compared  with  the  detected  arrival 
time  of  each  pulse.  If  a  detected  arrival  time  corresponds  to  one  of  the  determined  main 
multipath  arrival  times,  it  is  transferred  to  the  new  matrix  on  the  respective  row  of  the 
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pulse  from  which  it  was  detected  and  the  column  to  which  multipath  it  is  associated.  If 
no  arrivals  are  detected,  a  zero  is  inserted  in  the  matrix  as  a  place  holder.  The  size  of 
matrix  B  can  be  different  from  the  size  of  A,  depending  on  how  much  the  impulse 
response  of  the  channel  varies  between  the  transmission  times  of  the  n  pulses. 


j  -  multipath  arrival  index 
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Figure  9.  Example  of  matrix  B  generated  by  the  PDF  estimator  and  arrival  sorting 
algorithm.  B  is  the  matrix  of  sorted  arrival  times.  The  number  of  rows  in  B 
remains  the  same  as  the  number  of  rows  in  A  (Figure  8)  but  the  number  of 
columns  in  B  can  change  depending  on  the  number  of  multipath  arrivals  identified 
by  the  algorithm.  Zeros  are  inserted  where  no  arrivals  can  be  correlated  to  an 
identified  multipath  within  a  received  pulse. 


Now  aligned,  each  common  ensemble  of  arrivals  (individual  columns  of  matrix  B 
in  Figure  9)  is  processed  by  a  linear-regression  algorithm.  As  mentioned  in  Section  A  of 
this  chapter,  the  estimated  time  delay  measured  for  each  chirp  for  a  common  multipath 
arrival  and  the  known  RPI  can  be  used  to  perfonn  a  linear-regression  analysis.  A  straight 
line  is  fit  to  the  data  in  each  of  the  columns  of  B.  This  analysis  results  in  two  values:  Fj 

is  the  slope  of  the  best  fit  line  and  E-  represents  the  goodness  of  fit  of  this  line,  better 
defined  as  the  linear  correlation  coefficient  [14].  As  before,  the  slope  F-  represents  the 

skew  between  the  synchronizing  and  reference  clocks.  The  linear  correlation  coefficient 
Ej  is  used  in  the  next  step  of  the  algorithm  to  select  the  most  stable  multipath  arrival. 

The  final  step  of  the  process  identifies  which  set  of  arrivals  is  the  most  stable  and 
reliable  within  the  series  of  chirp  responses  analyzed.  This  decision  is  based  on  the 
calculated  value  of  E-  and  the  number  of  arrivals  detected  in  a  particular  set.  The  linear 

correlation  coefficient  is  calculated  based  on  the  standard  deviations  of  the  measured  and 
expected  time  delays  and  their  associated  mutual  covariance  [14].  The  value  of  the 
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mutual  coherence  can  vary  between  0  and  1,  where  1  represents  the  best  possible  fit  (i.e., 
all  the  points  are  on  the  line)  and  0  the  worst  possible  fit  (no  correlation  among  the 
arrivals  of  the  ensemble). 

Time  variability  within  a  specific  multipath  causes  the  arrivals  to  vary  in  the  time 
domain,  resulting  in  a  greater  deviation  of  the  arrivals  from  the  best-of-fit  line.  We  sort 
all  of  the  correlation  coefficients  Ej  j  =  l,2,...,m  and  find  the  largest  one.  This  most 

stable  multipath  arrival  is  chosen  to  be  the  one  with  the  highest  correlation  coefficient 
(i.e.,  closest  to  1),  and  this  arrival  has  the  least  deviation  from  its  best-of-fit  line. 

Another  metric  required  to  confirm  this  most  stable  multipath  arrival  is  the 
number  of  arrivals  included  in  the  ensemble.  A  low  number  of  arrivals  within  a 
multipath  arrival  ensemble  is  a  good  indication  by  itself  that  it  is  not  a  reliable  multipath 
arrival  (if  a  column  of  B  contains  several  zeros).  As  well,  it  can  negatively  influence  the 
result  of  the  linear  regression.  As  an  example,  perfonning  linear  regression  with  only 
two  arrivals  produces  a  perfect  correlation  coefficient  of  1 .  The  most  stable  multipath  is 
identified  by  the  number  of  multipath  detections  and  the  associated  time-variability  of 
that  multipath.  The  user  must  make  this  judgment. 

C.  STAGE  2:  OFFSET  CORRECTION 

Our  approach  to  correct  the  offset  bA  (as  defined  in  Equation  1)  of  the  skew- 
corrected  clock  compared  to  a  reference  is  different  than  THSL’s  as  it  does  not  assume 
that  the  propagation  between  the  two  nodes  is  the  same  as  shown  in  Figure  3.  The  reason 
channel  reciprocity  is  not  a  valid  assumption  in  the  underwater  environment  is  explained 
in  Chapters  II  and  III.  Such  a  situation  is  illustrated  in  Figure  10  by  separately  labeling 
the  one-way  propagation  delays. 
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Figure  10.  Two-way  message  exchange  of  timing  messages  between  the  synchronizing 
(A)  and  reference  (B)  nodes  where  the  individual  one-way  propagation  delay 
differs.  In  this  example  pA3  >  PH  A  .  After  [1]. 

The  definition  of  P  in  Equation  2  is  modified  as 

Pa,b+Pba  =  P •  (6) 

As  well.  Equations  3  and  4  are  modified  to  account  for  the  different  propagation 

delays 

T2~T\-PAfi=bA  (7) 

P4~P3~  Pb,A  ~  bA  ■  (^) 

The  basis  of  our  approach  is  the  two-way  message  exchange  used  by  most  clock- 
synchronization  protocols.  Its  distinction  lies  in  its  ability  to  differentiate  the  propagation 
delays  for  both  one-way  propagation  channels  instead  of  using  an  average  based  on  the 
round-trip  propagation  time.  A  chirp  is  sent  ahead  of  each  message  allowing  the  chirp 
response  of  the  one-way  communication  channel  to  be  detennined  by  the  receiver.  Both 
estimated  impulse  responses  are  correlated  at  the  synchronizing  node  A  to  detennine  the 
difference  in  propagation  delay  A p  between  the  two  of  them.  This  delay  Ap  is  the 
difference  between  the  time  corresponding  to  the  largest  peak  and  zero  lag.  This  method 
can  be  used  because  the  modem  always  synchronizes  on  the  largest  amplitude  arrival. 
The  offset  of  the  clock  can  be  calculated  using  this  propagation  delay  difference  and  the 
four  time-stamps  generated  by  the  two-way  message  exchange. 
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An  overall  system  diagram  of  this  adapted  two-way  message  exchange  is 
presented  in  Figure  1 1  and  is  further  explained  below. 


Figure  1 1 .  Adapted  two-way  message  exchange  for  calculating  the  offset  of  a  skew- 

corrected  local  clock  compared  to  a  reference. 

The  synchronizing  node  A  starts  the  offset  correction  process  by  sending  a  new 
pulse  ahead  of  the  offset  calculation  request  message  to  reference  node  B.  The  time  of 
the  start  of  the  transmission  is  time-stamped  Tx  and  recorded  by  node  A. 

As  it  propagates  through  the  channel,  the  pulse  is  transformed  as  explained  in 
Chapter  III.  The  yet-to-be-detennined  propagation  delay  is  denoted  as  pA  B  in  Figure  11. 
Time-stamp  T2  is  recorded  by  the  receiver  upon  first  detecting  the  arrival  of  the  pulse 
sent  by  node  A.  From  the  pulse  response,  the  instantaneous  one-way  impulse  response 
from  node  A  to  B  is  determined  using  a  matched  filter  and  identifying  the  peaks  of  the 
resulting  time-series.  This  impulse  response  is  required  later  in  the  algorithm  at  the 
synchronizing  node.  To  reduce  the  bandwidth  required  to  send  it  across  the  channel  for 
processing,  the  estimated  impulse  response  is  compressed  by  retaining  its  main 
characteristics,  which  are  the  amplitudes  of  the  largest  peaks  and  their  associated  time 
indices.  These  characteristics  are  represented  as  HAB  in  Figure  11. 

33 


The  reference  node  then  sends  the  same  new  reference  pulse  before  transmitting 
the  second  message  of  the  exchange.  Time-stamp  T3  is  recorded  concurrently  to  the 
beginning  of  the  transmission  of  the  second  message.  It  carries  the  two  time-stamps 
generated  by  node  B,  T2  and  T3 ,  as  well  as  the  compressed  impulse  response  information 
Hab  recorded  from  the  probe  signal  preceding  the  first  message. 

The  second  message  propagates  through  the  channel  from  node  B  to  A.  The 
message  is  transfonned  by  the  channel  as  discussed  in  Chapter  III.  The  yet  to  be 
detennined  propagation  delay  is  denoted  as  pBA  in  Figure  11.  Time-stamp  T4  is 

generated  by  the  receiver  upon  detecting  the  arrival  of  the  chirp  sent  by  node  B.  From 
this  received  chirp,  the  instantaneous  one-way  impulse  response  from  node  B  to  A 
(Hba)  is  estimated  and  its  infonnation  compressed  using  the  same  process  explained  in 
step  3  of  this  algorithm.  From  their  respective  recorded  characteristics,  the  main  arrivals 
of  both  estimated  impulse  responses  are  reconstructed  in  the  time-domain  from  HA3  and 
Hb  a  .  Their  cross-correlation  allows  for  the  time-delay  estimation  between  the  two 

signals.  Based  on  the  assumption  that  the  modem  will  synchronize  with  the  strongest 
arrival,  the  difference  in  propagation  delay  A p  (Equation  9)  will  be  the  time  index 
difference  between  the  highest  peak  of  the  cross-correlation  and  its  zero  lag: 

Pa,b-Pb,a  =  AP  =  z*-  (9) 

The  total  propagation  delay  P  of  both  messages  and  the  time-delay  difference 
A p  can  be  combined  to  evaluate  one  of  the  two  propagation  delays  ( pH  A  or  pAB). 
Adding  Equations  6  and  9  together,  we  get 

P  +  Ap  =  2pAB.  (10) 

The  one-way  propagation  delay  pAB  can  be  determined  from  Equation  10: 


P  +  Ap 

—  Pab  ■ 


(11) 


The  second  propagation  delay  (in  this  case  pH  A  )  can  be  evaluated  from 
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P  Pa,b  Pb,a  • 


(12) 


Either  one  of  the  one-way  propagation  delay  can  be  used  in  Equation  1 1  as  A p 
can  be  positive  or  negative.  Whichever  one-way  propagation  is  evaluated  in  Equation  11, 
the  remaining  one  can  be  detennined  using  Equation  12. 

The  last  step  of  the  algorithm  consists  of  evaluating  the  offset  bA  between  the  two 
nodes’  clock  using  Equations  7  or  8.  The  result  of  both  equations  should  be  the  same. 

D.  TIME  REVERSAL  APPROACH 

A  time  reversal  approach  was  investigated  as  a  way  to  evaluate  the  individual 
one-way  propagation  delays  of  the  two-way  message  exchanged  methodology.  This 
approach  initially  appeared  very  promising  since  one  of  its  benefits  is  the  elimination  of 
the  multipath  effect  of  the  communication  medium.  However,  one  of  the  assumptions  of 
time-reversal  is  the  requirement  of  reciprocity  within  the  channel  [16].  The  basis  of  our 
algorithm  is  the  existence  of  different  one-way  propagation  delays  resulting  in  a  non¬ 
reciprocal  channel.  Therefore,  a  time  reversal  approach  was  rejected  for  our  problem. 
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V.  COMPUTER  SIMULATIONS  STUDIES 


We  implemented  the  approach  described  in  Chapter  IV  in  Matlab.  The  Bellhop 
and  VirTEX  underwater  acoustic  propagation  models  (both  compatible  with  Matlab)  are 
used  to  simulate  a  simple  underwater  acoustic  communication  channel.  The  aim  of  the 
simulation  is  to  prove  our  conceptual  approach  prior  to  conducting  an  in-water 
experiment.  The  environmental  models,  the  simulated  environment’s  parameters,  the 
step-wise  results,  and  the  clock-correction  stages  (skew  and  offset)  are  explained  in  the 
following  sections. 

A.  BELLHOP 

Bellhop  [17]  is  a  ray-tracing  program  developed  by  Michael  B.  Porter  and 
distributed  as  part  of  the  Acoustic  Toolbox  software  package.  Based  on  ray  theory,  it  can 
produce  two-dimensional  acoustic  ray  tracing  for  a  given  channel  geometry  and 
underwater  environment.  Various  options  can  be  selected  and  inputs  given  to  the  model 
to  better  define  the  environment  and  provide  more  realistic  simulation.  Once  the  channel 
geometry  and  environment  are  defined,  Bellhop  provides  various  graphical  and  numerical 
outputs. 

B.  VIRTEX 

The  VirTEX  family  of  algorithms  was  designed  by  Peterson  and  Porter  [18].  The 
algorithms  build  on  the  results  obtained  by  a  Bellhop  simulation  and  add  the  effects  of 
multipath  and  Doppler  introduced  by  environmental  motion.  The  algorithm  computes  the 
time  series  that  would  be  observed  at  a  hypothetical  receiver  location. 

C.  DESCRIPTION  OF  THE  SIMULATED  ENVIRONMENT 

Our  simulation  environment  is  very  simple.  The  aim  of  our  simulation  is  not  to 
model  a  specific  environment  but  to  test  the  theory  of  our  algorithm.  The  environment, 
illustrated  in  Figure  12,  consists  of  a  transmitter  and  a  receiver,  both  situated  at  a  depth  of 
100  m  and  separated  by  a  distance  of  1  km.  The  sound-speed  profile  is  isothermal  for 
simplicity.  Time  variability  is  introduced  in  the  environment  by  randomly  changing  two 
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parameters:  the  orientation  of  the  transmitter  and  the  phase  of  a  simulated  wave  on  the 
sea  surface.  As  explained  in  Chapter  III,  the  orientation  of  the  beam  pattern  can  change 
the  path  which  has  the  most  energy.  As  well,  surface  waves  can  alter  the  delay-spread 
characteristics  of  the  impulse  response.  The  surface  wave  is  modeled  as  a  sinusoidal 
wave  with  phase  randomly  changed  throughout  the  simulation  to  provide  a  non-stationary 
surface  boundary. 


Normally,  both  nodes  of  the  system  will  operate  with  the  same  sampling 
frequency.  To  simulate  a  skew  difference  between  the  two  clocks,  the  sampling 
frequency  of  the  synchronizing  node  is  slightly  altered  when  recording  the  incoming 
chirp  time-series.  This  artificial  method  creates  a  known  skew  between  the  two  nodes’ 
clocks.  For  this  simulation,  the  sampling  frequency  at  the  receiver  is  reduced  from 
48,000  samples/sec  to  47,885  samples/sec.  This  reduction  in  sampling  frequency  results 
in  a  skew  of  417  ppm. 
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Figure  12.  Characteristics  of  the  simulated  environment.  From  the  top  left  corner:  the 
sound-speed  profile;  the  geometry  of  the  environment  and  its  associated  ray  paths; 
the  impulse  response  of  the  communication  channel. 
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D.  STEP-BY-STEP  SIMULATION  RESULTS 

This  section  presents  the  formal  results  of  the  simulations.  Figures  illustrating  the 
different  stages  of  the  algorithm  are  also  presented  to  help  visualize  some  of  the 
explanations  from  Chapter  IV. 

1.  Phase  1:  Skew  Correction 

A  5 -ms  LFM  chirp  based  on  a  5 -kHz  bandwidth  is  created  and  designated  as  the 
reference  chirp.  The  chirp  is  sent  through  the  channel  20  times  at  a  RPI  of  one  per 
second.  The  orientation  of  the  transmitter  and  phase  of  the  surface  wave  is  randomly 
changed  prior  to  each  transmission  and  remains  the  same  during  the  transmission  of  a 
single  chirp.  The  simulation  is  run  with  and  without  a  skew  difference. 

The  receiver,  knowing  the  sampling  frequency,  RPI  and  number  of  pulses  sent, 
separates  the  pulses  into  a  matrix  of  20  chirps  as  illustrated  in  Figure  13. 
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Chirp  responses 


Figure  13.  Ensemble  of  20  chirp  responses  with  no  clock  skew. 


Without  any  signal  processing,  time  variability  is  already  noticeable  in  Figure  13. 
Let  n  denote  the  number  of  probing  pulses  (chirps)  used  in  any  given  experiment,  n  =  20 
is  then  the  number  of  measurement  waveforms  displayed  in  Figure  13.  We  denote  these 

n  waveforms  by  {  V,(0}"=l  • 

We  compute  the  cross  correlation  R  (t)  between  the  reference  chirp  x(t)  and 
each  of  the  measurement  signals  y,(t)  to  create  the  ensemble  of  cross  correlations 

■[/?„,  (f)|  i  plotted  in  Figure  14. 
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Cross  correlated  responses  from  chirps 


Figure  14.  Ensemble  of  20  cross-correlated  received  chirp  responses  with  no  clock  skew 

present. 

Note  that  the  signals  in  Figures  13  and  14  contain  m  =  5  multipath  arrivals  that 
make  up  the  chirp  response.  Let  j  denote  the  “arrival  index,”  or  the  index  on  the  m 

r  i  m 

multipath  arrivals  that  occur  in  a  given  single  measurement 

waveform,  j  =  1, 2, . . . ,  m  . 

The  peaks  displayed  in  Figure  15  are  estimated  from  the  cross-correlation  signals 
displayed  in  Figure  14.  First,  the  individual  arrivals  in  the  cross-correlation  waveforms 
are  selected  to  be  the  multipath  arrivals  that  exceed  a  small  threshold  chosen  by  the  user. 
For  each  of  the  arrivals  in  a  given  correlation  wavefonn  Rxy  (r) ,  we  find  the  time  T* 

corresponding  to  the  time  at  which  a  peak  occurs  by  choosing  the  time  lag  t  such  that 
(^ * )  =  max  { (T )}  . 
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Figure  15.  The  peak  times  for  an  experiment  in  which  clock  skew  is  present:  The  peak 
times  Uji  derived  from  the  ensemble  of  n  =  20  measurement  signals  depicted  in 

Figure  14  are  plotted  in  this  figure.  Note  that  each  of  the  m  -  5  arrival  wavelets 

a.j ,  j  =  1, 2, . . .,  m  has  a  corresponding  ensemble  £/;  =  [n;!  J  ,  j  =  1, 2, . . . ,  m 

associated  with  it. 

We  see  in  Figure  14  the  ensemble  of  n  waveforms  and  in  Figure  15  the  ensemble 
of  peaks  of  waveforms  that  result  from  the  n  =  20  probing  pulses  from  the  transmitter. 
After  peak  detection,  we  see  in  Figures  15  and  16  that  each  multipath  arrival  has  a 
corresponding  ensemble  of  n  =  20  peaks  associated  with  it.  The  peak  arrival  times  shift 
over  time  due  to  the  skew  between  the  transmitter  and  receiver  clocks.  An  example  of  an 
ensemble  of  waveforms  with  zero  skew  is  depicted  in  Figure  16. 
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Figure  16.  The  peak  times  for  an  experiment  in  which  zero  clock  skew  is  present:  The 
peak  times  u  fi  derived  from  the  ensemble  of  n  =  20  measurement  signals  depicted 

in  Figure  14  are  plotted  in  this  figure.  Note  that  each  of  the  m  =  5  arrival 
wavelets  a .,  j  =  1, 2, . . m  has  a  corresponding  ensemble 

Uj  =  [n>( }  ,j  =  l,2,...,m  associated  with  it. 


Let  i  denote  the  peak  index  corresponding  to  each  of  the  n  peaks  that  occur  for  a 
single  multipath  arrival  in  a  given  ensemble  constructed  from  the  data  displayed  in 
Figures  15  and  16.  Because  n  probe  pulses  create  an  ensemble  of  n  waveforms,  there 
are  n  peaks  in  an  ensemble  corresponding  to  the  y*  multipath  arrival  in  Figures  15  and 

16,  i  =  l,2,...,n  .  We  denote  the  zth  peak  time  for  the  ylh  multipath  arrival  by  «/7 .  We 

denote  the  ensemble  of  peak  times  for  they*  multipath  arrival  by  Uj  =  {uji} 
j  =  1,2,..., m  . 

For  each  multipath  arrival  a  .,  j  =  1, 2,...,  m  ,  we  estimate  a  PDF  of  peak  arrival 
times  using  the  ensemble  Uj  =  |n;,  j  ,  j  =  1, 2, . . . ,  m  .  We  use  a  kernel  density  estimator 


42 


with  a  Gaussian  kernel  to  build  the  PDF  estimates  [19].  Once  all  m  PDFs  are  estimated 
individually,  we  plot  them  on  a  single  graph  as  depicted  in  Figure  17.  We  see  that  in 
Figures  15,  16  and  17  there  are  m  —  5  multipath  arrivals. 

Individual  PDFs  of  the  arrivals  for  each  path  (blue)  and  impulse  response  of  the  channel  (red) 


Figure  17.  The  scaled  simulated  impulse  response  (a  set  of  Kronecker  delta  functions)  of 
the  channel  is  plotted  in  red.  The  estimated  PDFs  of  the  individual  arrival 
wavelets  are  plotted  on  a  single  graph  and  appear  in  blue.  Note  that  the  delta 
functions  occur  at  the  times  at  which  the  chirp  pulses  begin,  so  the  centers  of  the 
PDF  estimates  are  shifted  slightly  to  the  right  of  delta  functions. 

The  shape  of  each  PDF,  as  illustrated  in  Figure  17,  is  representative  of  the 
stability  of  each  signal  propagation  path.  The  direct  paths  and  bottom-bounce  paths  are 
not  affected  by  the  time  variability  of  the  environment,  which  is  only  present  at  the  water 
surface  interface.  These  unaffected  paths  are  characterized  by  sharp  and  narrow  PDFs 
whereas  the  surface-bounce  paths  are  wide  and  low.  To  demonstrate  the  time-variability 
of  the  environment,  the  non-Doppler  affected  estimated  impulse  response,  as  shown  in 
Figure  12,  is  superimposed  on  Figure  17.  It  is  important  to  note  that  the  amplitudes  of 


43 


the  PDFs  are  not  related  to  the  amount  of  energy  propagated  through  the  associated  path. 
The  third  arrival  shown  in  Figure  17  is  a  good  example  of  the  case  in  which  the  weakest 
signal  path  is  very  stable. 


Individual  PDFs  of  the  arrivals  for  each  path  (blue)  and  impulse  response  of  the  channel  (red) 
x  io-3  zero  clock  skew  present 


Figure  18.  The  scaled  simulated  impulse  response  (a  set  of  Kronecker  delta  functions)  of 
the  channel  is  plotted  in  red.  The  estimated  PDFs  of  the  individual  arrival 
wavelets  are  plotted  on  a  single  graph  and  appear  in  blue.  Note  that  the  delta 
functions  now  occur  in  the  middle  of  the  chirp  pulses,  as  zero  skew  is  present. 

Figure  18  is  a  replica  of  Figure  17  but  with  zero  skew  present.  Since  the  1st  and 
3rd  arrivals  are  not  affected  by  the  time  variability  of  the  channel  and  zero  skew  is 
present,  their  respective  PDFs  are  now  represented  by  Kronecker  delta  functions.  In 
contrast  to  Figure  17,  the  delta  functions  now  occur  in  the  middle  of  the  chirp  pulses. 

A  least-squares  linear  regression  algorithm  is  applied  to  each  ensemble  of  arrivals. 
The  result  of  this  linear  regression  is  the  clock  skew  difference  F- ,  from  which  the 

coefficient  of  determination  E,  can  be  calculated  as  explained  in  Chapter  IV.  Since  the 

tolerance  of  a  clock  is  usually  measured  in  ppm,  Equation  13  can  be  used  to  convert  the 
skew  to  Fpj  to  compare  the  measured  value  with  the  published  tolerance  of  a  clock: 
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\,000,000(Fr\)=Fpr 


(13) 


The  simulation  results  for  the  first  phase  of  the  clock-synchronization  algorithm 
are  illustrated  in  Table  1.  The  algorithm  chooses  the  first  path  as  the  most  reliable 
(shaded  arrival  in  Table  1)  as  it  has  the  highest  correlation  coefficient  (perfect  in  this 
case).  This  is  expected,  since  both  the  first  and  third  arrivals  do  not  interact  with  the  sea 
surface  and,  therefore,  no  time  variability  is  introduced  in  their  associated  signal 
propagation  path  (direct  and  single  bottom  bounce)  as  illustrated  in  Figure  12. 


Table  1.  Skew  correction  simulation  results.  The  shaded  result  corresponds  to  the  most 
stable  arrival  set  as  selected  by  the  algorithm  (the  set  with  the  largest  correlation 

coefficient). 


Simulation 
and  set  clock 
skew 

Arrival 

Ensemble 

Wj) 

Skew  ( Fj ) 

Correlation  coefficient 
(■ EJ ) 

Skew  in 
ppm 

( FPj ) 

#1 

No  clock 
skew 
(0  ppm) 

1 

1.000000000000000 

1.000000000000000 

0 

2 

1.000087954260652 

0.999999926902282 

88 

3 

1.000000000000000 

1.000000000000000 

0 

4 

1.000289301378446 

0.999999011969708 

289 

5 

1.000683255012531 

0.999995976887908 

683 

#2 

With  clock 
skew 

(417  ppm) 

1 

1.000416666666667 

1.000000000000000 

417 

2 

1.000504620927318 

0.999999926963154 

504 

3 

1.000416666666667 

1.000000000000000 

417 

4 

1.000705968045113 

0.999999012792314 

705 

5 

1.001099921679198 

0.999995980236109 

1099 

After  artificially  introducing  a  clock  skew  in  the  simulation,  the  algorithm 


successfully  detennines  the  skew  as  417  ppm,  which  is  the  value  expected  from  the 
reduction  in  sampling  frequency  at  the  synchronizing  node  (48,000  samples/sec  to  47,885 
samples/sec). 

This  simulation  demonstrates  that  time  variability  can  affect  the  measurement  of 
the  skew  between  two  clocks.  The  skew  values  calculated  from  the  ensembles  of  the 
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multipath  arrival  times  affected  by  time  variability  varies  greatly.  In  the  worst  case 
(ensemble  #5)  the  skew  value  in  ppm  is  nearly  three  times  the  real  skew  value. 

2.  Stage  2:  Offset  Correction 

The  same  environment  and  reference  chirp  as  in  stage  1  are  used  to  simulate  our 
approach  to  the  clock-offset  correction.  The  simulation  is  designed  to  follow  the  steps 
explained  in  Chapter  IV  and  illustrated  in  Figure  11.  It  focuses  on  the  calculation  of  the 
difference  A p  in  propagation  time  between  the  two  messages  sent  in  the  two-way  timing 
message  scheme.  Once  this  parameter  is  determined,  the  other  unknown  values 
( P,pAfi,pBA  and  bA)  can  be  found  using  the  four  recorded  time-stamps  and  Equations  2, 

1 1  and  12.  An  example  of  how  these  parameters  can  be  calculated  with  these  equations 
is  illustrated  in  Figure  19. 

Synchronizing  Reference 

Node  A  Node  B 

Ti=  50 

r4  =  73 

|  Ap=-2  | 

Figure  19.  Example  of  time-stamp  exchange  to  detennine  the  clock  offset  bA  between 

nodes  A  and  B. 

From  Equation  2,  the  total  propagation  time  P  is  calculated  as 

(T4  -T1)-(T3-T2)  =  (63-  50) -(40  -  35)  =  18  =  P.  (14) 

Then,  using  Equation  1 1 ,  and  knowing  P  and  A p  ,  we  get 
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(15) 


P  -  Ap 
2 


18 -(-2) 
2 


=  10  = 


Pas- 


The  last  propagation  delay  pB  A  can  be  calculated  from 

Pab  =  18  - 10  =  8  =  pBA.  (16) 

Finally  the  offset  bA  can  be  detennined  using  either  Equations  7  or  8. 

To  test  our  approach,  the  simulation  must  create  an  environment  where  a  Ap  will 
exist  between  the  two  message  propagation  delays.  Assuming  that  the  modem  (node) 
picks  the  strongest  arrival  to  synchronize  with,  we  see  that  the  simulation  must  vary  the 
location  of  this  arrival  within  the  impulse  response  of  the  two  message  transmissions  to 
create  a  Ap .  The  parameters  used  in  Figure  4,  where  the  transmitter  was  tilted  from  0  to 
30  degrees,  were  used  again  for  this  simulation  since  they  illustrate  the  desired  temporal 
displacement  of  the  strongest  arrival.  The  steps  of  the  simulation,  illustrating  the  various 
signal  processing  steps,  are  described  in  this  section. 

The  synchronizing  node  A  initializes  the  clock  offset  correction  phase  following 
the  determination  of  the  clock  skew  in  stage  1.  It  sends  a  reference  chirp  and  an  offset 
correction  request  message  to  the  reference  node  B  while  also  creating  and  storing  a  time 
stamp  7,  when  the  transmission  is  initiated. 

The  received  chirp  response,  illustrated  in  Figure  20,  is  processed,  where  its  main 
characteristics  (peak  amplitudes  and  corresponding  time  index)  are  extracted  in  order  to 
limit  the  amount  of  information  transmitted  across  the  channel.  Time  stamps  T2  and  T3 
are  respectively  created  by  the  reference  node  when  it  receives  the  first  message  and 
transmits  the  second  one.  The  difference  between  these  two  time  stamps  represents  the 
processing  or  handling  of  the  request  by  the  reference  node.  The  reference  node 
compiles  the  two  time-stamps  and  the  characteristics  of  the  analyzed  estimated  impulse 
response  into  a  message,  which  is  sent  back  to  the  synchronizing  node  following  the 
transmission  of  another  chirp. 
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Figure  20.  LFM  chirp  response  received  by  the  reference  node  B. 

The  synchronizing  node  receives  the  chirp  response  and  processes  it  in  the  same 
manner  as  did  the  reference  node.  This  results  in  a  basic  estimated  impulse  response  for 
each  one-way  transmission  between  the  two  nodes,  illustrated  in  Figure  21. 
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Figure  21. 


Simplified  estimated  impulse  responses  from  the  two  transmissions 


These  two  estimated  impulse  responses  are  cross-correlated  resulting  in  the  cross 
correlation  signal  plotted  in  Figure  22. 
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Figure  22.  Cross-correlation  of  both  simplified  impulse  responses.  The  red  line  indicates 
the  middle  of  the  cross-correlation.  The  distance  between  the  highest  peak  and  the 
middle  of  the  cross-correlation  represents  A p  . 

The  delay  between  the  highest  peak  and  the  center  of  the  cross-correlation  (time 
index  0)  is  the  difference  in  propagation  delay  A p  between  the  two  messages  sent.  The 
measured  time  series  simulated  by  VirTEX  is  referenced  in  the  time  domain  to  the 
beginning  of  the  transmission,  as  shown  in  Figure  20,  which  allows  the  determination  of 
the  “true”  propagation  delays  pB  A  ,  pAB  and  A p  using  Equation  9. 
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Table  2.  True  and  calculated  values  of  the  difference  in  propagation  delays  between  the 

two  messages  during  the  two-way  exchange. 


True  PAS 

0.6664  sec 

True  pBA 

0.6536  sec 

True  Ap 

-0.0128  sec 

Measured  Ap  from  cross-correlation 

-0.0128  sec 

As  stated  at  the  beginning  of  this  section,  once  Ap  is  detennined,  the  clock  offset 
is  easily  calculated.  The  simulation  proves  the  theory  of  our  approach  since  it 
successfully  and  precisely  determines  the  difference  in  propagation  delays  and  clock 


offset. 
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VI.  DEL  MONTE  LAKE  EXPERIMENT 


A.  OBJECTIVE 

The  aim  of  the  experiment  is  to  test  and  prove  the  algorithms  described  in  Chapter 
IV  in  an  actual,  time-variant,  underwater  environment.  Only  the  first  stage  (skew 
correction)  is  tested  during  the  experiment.  The  second  stage  requires  resources  not 
available  to  this  project,  such  as  access  to  the  firmware  of  the  Teledyne  Benthos  acoustic 
modems. 

Del  Monte  Lake,  situated  on  the  campus  of  the  Naval  Postgraduate  School  is 
selected  for  its  accessibility  and  minimal  environmental  studies  requirements  and 
approvals.  The  general  concept  of  the  experiment  consists  of  transmitting  and  recording 
a  series  of  chirps  between  two  acoustic  modems.  The  recorded  data  are  downloaded 
from  an  on-board  digital  recorder  and  processed  separately.  During  the  post-experiment 
processing,  the  arrivals  of  chirp  responses  are  tracked  and  processed  to  determine  the 
skew  difference  between  the  clocks  of  the  two  acoustic  modems  (nodes).  The  lake 
environment  on  the  day  of  the  experiment,  the  equipment  used,  the  test  procedures, 
selection  of  parameters  and  an  evaluation  of  the  results  are  explained  in  more  detail  in  the 
following  section. 

B.  ENVIRONMENTAL  DESCRIPTION 

The  experiment  was  conducted  during  the  afternoon  of  2  May  2012  on  the  Del 
Monte  lake.  A  satellite  picture  [20]  of  the  area  is  illustrated  in  Figure  23. 
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Figure  23.  Satellite  picture  of  Del  Monte  Lake  and  locations  of  the  acoustic  modems 

during  the  experiment  (from  [20]). 

It  was  a  warm  sunny  day  with  reported  winds  of  24  km/h  [21].  The  wind  was 
strong  enough  to  generate  waves  on  the  surface  of  the  lake  and  made  it  difficult  to 
maintain  a  stationary  location  with  the  small  aluminium  boat.  From  a  previous 
experiment  on  the  lake,  the  depth  of  the  water  was  found  to  be  approximately  2  m  over 
most  of  the  lake.  The  bottom  of  the  lake  consists  of  thick  mud.  The  water  temperature 
was  not  measured  but  was  assumed  to  be  isothermal  through  the  water  column  due  to  the 
very  shallow  depth.  The  water  fountain  (not  shown  on  Figure  23)  was  turned  off  to 
reduce  the  amount  of  noise  present  in  the  lake. 

C.  EQUIPMENT  DESCRIPTION 

1.  Aluminum  Boat  and  Motor 

A  13  foot  aluminum  boat  equipped  with  a  Minn  Kota  ENDURA  Pro  C2  transom- 
mounted  trolling  freshwater  electric  motor  was  used  to  conduct  the  experiment.  It  was 
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used  to  deploy  the  acoustic  modem  and  served  as  a  working  platform  on  which  the  deck 
box  and  laptop  were  operated  to  perform  the  experiment. 

2.  Teledyne  Benthos  SM-75  SMART  Modem 

A  Teledyne  Benthos  SM-75  Smart  Modem  Acoustic  Release  Technology 
(SMART)  [5]  was  refurbished  and  fitted  with  an  onboard  digital  recorder,  as  illustrated  in 
Figure  24. 


Figure  24.  Teledyne  Benthos  SM-75  SMART  modem  (from  [5]). 

The  modem  operates  in  the  low-frequency  (LF)  range  of  9  -14  kHz  and  uses  a 
sampling  frequency  of  10240  samples/sec.  The  digital  recorder  is  equipped  with  a  4-GB 
Secure  Digital  High  Capacity  (SDHC)  card,  which  records  the  baseband  received  signal 
by  the  modem  in  minute-long  .wav  or  .pern  files.  The  gain  of  the  received  signal  is 
continually  adjusted,  using  an  automatic  gain  controller  (AGC),  to  keep  the  signal  level 
constant  in  the  recording.  The  AGC  information  is  also  recorded  as  part  of  the  baseband 
signal  by  the  digital  recorder.  A  command  is  sent  to  the  SM-75  to  start  and  stop  the 
digital  recording,  which  allows  the  recording  of  selected  desired  events.  Once  recovered, 
the  SM-75  can  be  dismantled  and  the  SDHC  card  containing  the  recorded  data  can  be 
removed.  These  recordings  are  then  uploaded  to  a  stand-alone  computer  for  post- 
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processing  and  analysis.  For  better  results,  Teledyne  Benthos  provides  a  Matlab  routine 
which  uses  the  AGC  infonnation  of  a  recording  to  restore  the  signal  to  its  original  level 
and  continuity. 

The  LF  beam  pattern  of  the  SM-75  is  illustrated  in  Figure  25.  This  beam  pattern 
is  used  in  Chapter  II  to  demonstrate  the  effect  of  a  transducer  rotation  on  the  impulse 
response  of  the  channel. 


10  0  350 


Figure  25.  SM-75  LF  (9-14  kHz)  beam  pattern  (from  [5]). 


During  the  experiment  the  SM-75  was  deployed  and  anchored  at  station  A,  as 
shown  in  Figure  23,  at  a  depth  of  approximately  1  m.  It  functioned  as  the  synchronizing 
node  B,  recording  the  series  of  chirp  response  transmitted  by  the  reference  node  A  (deck 
box). 


3.  UDB-9400  Universal  Deck  Box 

As  illustrated  in  Figure  26,  the  deck  box  is  an  acoustic  modem  where  the 

transducer  and  modem  electronics  have  been  separated.  Several  interfaces,  including  a 
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touch  screen  display,  are  combined  with  the  modem  in  a  portable  case,  which  allows  a 
user  to  perform  network  administration  within  the  network.  The  transducer  is  linked  to 
the  portable  case  via  a  long  cable  and  can  be  deployed  over  the  side  of  a  boat  or  platform. 


Figure  26.  UDB-9400  Universal  Deck  Box,  including  its  deployable  transducer  (left) 

(from  [5]). 

The  modem  part  of  the  deck  box  uses  the  same  communication  parameters  used 
by  the  SM-75  as  described  in  the  previous  section.  The  deck  box  is  not  equipped  with  a 
digital  recorder  but  still  has  a  limited  flash  memory  module,  where  files  can  be  uploaded. 
For  this  experiment,  .wav  audio  files  of  the  various  series  of  chirps  were  uploaded  to  the 
memory  from  which  they  could  be  played.  Instead  of  using  the  touch  screen  to  input 
modem  commands,  the  deck  box  was  connected  via  serial  port  (RS-232)  to  a  laptop 
running  the  PROCOMM  software.  PROCOMM  provides  an  easier  user  interface  and 
allows  the  creation  of  script  files  that  made  the  conduct  of  the  experiment  easier  and 
faster  to  complete.  The  same  laptop  can  also  be  interfaced  to  the  SM-75  or  acoustic 
modem. 

The  deck  box  was  used  as  the  reference  node  A  during  the  experiment.  Various 
series  of  chirps  were  transmitted  at  a  distance  of  2  m  (station  B  on  Figure  23)  and  50  m 
(station  C  on  Figure  23).  The  transducer  was  attached  half-way  down  a  3  m  long  metal 
pole  which  was  lowered  and  embedded  in  the  mud  at  the  bottom  of  the  lake.  The 
transducer  was  kept  stable  this  way  in  the  middle  of  the  water  depth  at  approximately  1 
m. 
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D.  EXPERIMENTAL  PARAMETERS 


The  parameters  of  the  experiment  had  to  be  chosen  to  accommodate  the 
characteristics  of  the  test  environment.  Underwater  acoustic  modems  are  not  normally 
deployed  in  a  shallow  fresh  water  lake  and  limited  information  was  available  about  any 
previous  experiments  in  such  environments.  To  mitigate  the  risk  associated  with  this 
unknown  environment,  various  parameters  were  tested  during  the  experiment.  These 
parameters  are  explained  in  the  next  sections. 

1.  Distance  Between  the  Nodes 

The  shallow  depth  of  the  lake  makes  the  communication  channel  prone  to 
multipath.  The  simultaneous  arrival  of  many  ray  paths  at  the  same  time  can  make  the 
retrieval  of  any  distinct  signal  impossible  due  to  temporal  superposition.  For  a  2  m  water 
depth,  the  maximum  propagation  delay  between  the  direct  and  surface  or  bottom  bounce 
paths  is  obtained  by  separating  the  nodes  by  2  m.  Assuming  an  isothermal  water  column 
and  a  constant  sound  of  speed  of  1530  m/s,  the  modeled  time-delay  difference  between 
the  two  main  paths  is  calculated  to  be  0.0005  s  compared  to  0.0001  s  when  the  nodes  are 
separated  by  50  m.  These  two  distances  were  used  during  the  experiment.  The  ray  paths 
and  impulse  response  for  nodes  separated,  respectively,  by  2  m  and  50  m  were  modeled 
using  Bellhop  and  are  shown  in  Figure  27. 
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Figure  27.  Bellhop  model  of  the  ray  paths  and  impulse  responses  in  Del  Monte  Lake  with 

2  m  and  50  m  node  separation. 


2.  Probe  Signal  Characteristics 

The  design  parameters  of  the  probe  signals  are  important  as  they  determine  how 
well  each  ray  path  can  be  distinguished  at  the  receiver.  It  is  generally  accepted  that  the 
LFM  chirp  is  a  good  choice  as  a  probe  signal  [10].  The  experiment  used  both  LFM 
chirps  and  CW  pulses  in  order  to  give  some  measure  of  the  advantage  of  using  the  LFM 
chirp.  The  frequency  of  the  chirp  swept  over  the  range  of  0  to  5000  Hz  which  was  later 
modulated  by  the  modem  to  its  operational  bandwidth  of  9  to  14  kHz. 

The  length  of  the  probe  signal  affects  the  amount  of  temporal  superposition  at  the 
receiver.  As  an  example,  when  the  nodes  are  separated  by  2  m,  the  time-delay  difference 
between  the  two  main  paths  is  such  that  a  0.5  ms  probe  signal  will  not  be  superimposed. 
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At  a  sampling  frequency  of  10240  Hz,  a  0.5  ms  probe  signal  corresponds  to  only  6 
samples  in  length.  Based  on  the  fact  that  the  ambiguity  function  of  the  chirp  is  inversely 
proportional  to  its  length  [22],  a  longer  chirp  should  provide  better  time  resolution. 
During  the  experiment,  chirps  of  0.5,  1,  2,  5,  10,  25  and  50  ms  duration  were  used  and 
compared  in  the  post-processing  analysis. 

The  pulse  repetition  interval  was  selected  and  remained  at  250  ms  throughout  the 
experiment.  This  value  was  arbitrarily  selected  as  it  seemed  short  enough  to  respect  the 
requirement  to  minimize  the  utilization  of  the  network  bandwidth  when  perfonning  a 
clock  synchronization  procedure  (as  explained  in  Chapter  II).  Another  factor  that  was 
considered  was  that  the  interval  be  short  enough  so  that  the  AGC  does  not  reset  in- 
between  each  probe  transmission. 

The  final  parameter  was  the  number  of  pulses  sent  in  a  given  series.  This  was 
selected  to  be  30  based  on  the  experimental  result  of  Syed  [1].  During  that  experiment,  it 
was  detennined  that  25  pulses  was  the  optimum  number  of  pulses  before  observing  a 
diminishing  benefit  in  increasing  the  number  of  pulses.  Since  Syed’s  experiment  was 
conducted  over  the  air,  our  number  of  pulses  was  increased  to  30  to  ensure  that  25  of 
them  would  be  received  in  the  less  favorable  underwater  environment. 

E.  TEST  PROCEDURES 

The  experiment  was  executed  in  three  stages,  which  translate  into  the  preparation, 
the  acquisition  of  measurements,  and  the  analysis.  The  procedures  associated  with  each 
stage  are  described  in  this  section. 

1.  Preparation 

In  preparation  for  the  experiment  the  above-stated  experimental  parameters  were 
chosen.  A  series  of  30  identical  pulses  was  generated  using  Matlab  and  then  converted  to 
a  .wav  file  format.  Each  file  was  uploaded  to  the  deck  box’s  internal  flash  memory  via 
the  serial  port  using  a  standard  laptop  computer.  Once  loaded  in  the  deck  box’s  memory, 
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the  .wav  file  could  be  played  through  the  modem  and  its  transducer.  A  PROCOMM 
script  file  was  written  to  execute  the  experiment  more  efficiently  and  with  fewer  manual 
errors. 

Before  proceeding  to  the  lake,  an  over-the-air  experiment  was  conducted  in  the 
Seaweb  Lab  to  confirm  that  the  equipment  worked  as  planned.  The  original  script  called 
for  the  digital  recorder  to  be  turned  on  and  off  between  the  transmission  of  each  series  of 
pulses,  which  sometimes  caused  the  SM-75  to  reboot  randomly.  The  script  was  modified 
so  that  the  digital  recorder  was  turned  on  at  the  beginning  of  the  experiment  and  only 
turned  off  once  every  series  of  pulses  was  transmitted.  This  solution  was  effective  and 
only  one  measurement  was  lost  due  to  a  SM-75  reboot. 

2.  Experimental  Measurements 

The  equipment  was  loaded  into  the  aluminum  boat  and  transported  to  the 
deployment  stations  as  shown  in  Figure  23.  The  SM-75  was  deployed  and  its  location 
recorded  using  a  handheld  GPS.  The  2-m  separation  between  the  SM-75  and  the  deck 
box  transducer  was  confirmed  by  the  handheld  GPS  and  by  the  ranging  routine  included 
in  the  onboard  software  of  the  Teledyne  Benthos  modems.  Every  pre-recorded  series  of 
probe  signals  was  transmitted  twice  at  both  the  2-m  (experiment  #1  and  #2)  and  50-m 
(experiment  #3  and  #4)  stations,  for  a  total  of  four  experiments.  No  additional 
experiments  were  perfonned  due  to  the  drained  battery  of  the  laptop  which  was  essential 
to  the  conduct  of  the  experiment.  The  SM-75  was  recovered  and  the  equipment  was 
returned  to  the  Seaweb  Lab.  The  recordings  from  the  SM-75  were  transferred  to  a  laptop 
for  analysis. 


3.  Post-experiment  Analysis 

A  log  file  was  created  along  with  each  recording  by  the  digital  recorder.  The  log 
files  identify  every  .wav  file  of  duration  of  one-minute  long  or  less  associated  with  it. 
Time-stamps  were  also  generated  by  the  modem’s  internal  clock  and  annotated  in  the  log 
marking  the  beginning  and  ending  of  each  files.  As  stated  earlier  in  this  chapter, 
Teledyne  Benthos  provides  a  Matlab  routine  that  reconstructs  the  recording  based  on  the 
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log  file.  With  this  routine,  each  recording  was  processed  and  the  original  non-AGC  time 
series  were  then  available  in  a  Matlab  format. 

Each  recording  contained  multiple  transmissions  of  pulse  series.  These  series 
were  easily  identifiable,  as  shown  in  Figure  28,  and  were  individually  extracted  to  be 
processed  separately. 


Figure  28.  AGC-corrected  recording  where  all  series  of  probe  signals  were  received 

when  the  nodes  were  separated  by  50  m  (experiment  #3). 

Each  individual  series  of  probe  signals  was  processed  separately  using  the  same 
algorithms  previously  described  in  Chapter  IV  and  V. 
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F.  RESULTS 

The  cumulative  results  of  the  analysis  are  illustrated  in  Table  3.  Following  Table 
3,  a  selection  of  these  results  is  expanded  to  demonstrate  some  of  the  noteworthy 
findings. 


Table  3.  Estimated  skew  (in  ppm)  and  its  related  error  for  each  different  series  of  probe 
signals  transmitted  during  each  experiment.  An  average  of  the  error  (in  ppm)  for 
each  different  series  of  probe  signals  is  also  included. 


Estimated  skew  and  error  in  ppm  for  each  experiment  (all  values  are  in  ppm) 

Probe  Signal 

Experiment 
#1  (2m) 

Experiment 
#2  (2m) 

Experiment 
#3  (50m) 

Experiment 
#4  (50m) 

Average 

Error 

0.5ms  LFM 

-42  ±  10 

-58  ±  15 

-33  ±  14 

-28  ±6 

±11.25 

1ms  LFM 

-42  ±  10 

-21  ±25 

4  ±  16 

-55  ±  14 

±16.25 

2ms  LFM 

-28  ±  12 

-53  ±  17 

-45  ±3 

-35  ±4 

±9 

5msLFM 

-38  ±  11 

-34  ±11 

-33  ±11 

-29  ±4 

±9.25 

lOmsLFM 

- 

-44  ±  50 

-20  ±  14 

-28  ±5 

±23 

25msLFM 

-41  ±  14 

-54  ±  22 

-31  ±6 

-51  ±  10 

±13 

50ms  LFM 

60  ±27 

-32  ±  27 

-34  ±11 

11  ±  15 

±20 

0.5ms  CW 

11  ±34 

5  ±46 

-56  ±25 

-48  ±3 

±27 

1ms  CW 

-46  ±8 

5  ±33 

-8  ±31 

-15  ±  15 

±21.75 

2ms  CW 

70  ±38 

-44  ±4 

-36  ±  20 

-40  ±  27 

±22.25 

5ms  CW 

-42  ±  33 

9  ±35 

-13  ±  16 

-59  ±21 

±26.25 

10ms  CW 

-30  ±32 

-39  ±39 

-28  ±4 

-16  ±26 

±25.25 

The  digital  recorder  rebooted  during  the  playing  of  the  10ms  LFM  series  of  pulses 


in  experiment  #1,  therefore,  no  results  are  shown  for  this  event  in  Table  3. 

An  initial  glance  seems  to  demonstrate  a  high  variance  among  the  results.  To 
better  analyze  this  variance,  the  uncertainty  of  each  skew  measurement  is  calculated  [14]. 
The  following  Figures  and  Tables  examine  the  source(s)  of  this  variance  and  provide  a 
higher  degree  of  confidence  in  the  experimental  results. 

From  Table  3,  the  2  ms  LFM  chirp  is  found  to  be  the  superior  probe  for  this 
environment  based  on  its  lower  average  error  of  ±9  ppm.  In  Figures  29  and  30, 
respectively,  the  first  1 5  matched- filtered  responses  of  the  2ms  chirp  when  tested  at  2  m 
(experiment  #2)  and  50  m  (experiment  #3)  are  shown.  These  figures  are  used  as  a 
comparison  with  other  chirps  that  were  identified  to  be  less  effective. 
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Figure  29. 


Figure  30. 


First  15,  2ms  chirps  (experiment  #2) 


Time  sample  index 


1200  15 


Chirp  index  (i) 


First  fifteen  2  ms  long  LFM  chirp  responses  at  2  m  (experiment  #2). 


First  15,  2ms  LFM  chirps  (experiment  #3) 
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Looking  at  Figures  29  and  30,  we  see  that  the  first  detected  and  most  likely  the 
direct  path  is  easily  identifiable  and  looks  fairly  stable  over  the  measurements.  The 
following  arrivals  are  not  as  easy  to  distinguish  but  a  few  peaks  are  estimated  to  be 
associated  with  a  multipath  arrival  by  the  algorithm  as  shown  in  Tables  4  and  5.  The  data 
presented  in  these  tables  are  the  time  index  where  a  peak  is  detected.  Perfonning  linear 
regression  on  each  of  these  paths  allows  the  determination  of  the  clock  skew  and  its 
correlation  coefficient,  which  is  used  to  identify  the  most  stable  path.  The  shaded  column 
highlights  the  most  stable  path  based  on  its  correlation  coefficient.  In  both  cases,  the  first 
path  has  the  highest  correlation  coefficient,  which  is  consistent  with  Figures  29  and  30. 
The  first  detected  path  clearly  shows  a  slow  shift  of  the  detected  peak,  indicating  the 
presence  of  skew  between  the  clocks  in  the  two  nodes.  Respectively  following  Tables  4 
and  5  are  Figures  31  and  32,  which  illustrate  the  same  data  presented  in  the  tables.  These 
figures  clearly  show  the  time  variability  of  the  multipath  arrivals. 


Table  4.  Detailed  results  of  the  2  ms  LFM  pulse  received  during  experiment  #2,  including 
the  correlation  coefficient,  clock  skew  (in  absolute  and  ppm  values)  and  its 

associated  error. 


Analysis  of  2  ms  LFM  pulse  received  during  experiment  #2  (2  m  separation) 

Pulse  Index 

Time  index  of  detected  peaks  for  each  multipath 

Detected 
path  #1 

Detected 
path  #2 

Detected 
path  #3 

Detected 
path  #4 

Detected 
path  #5 

1 

362 

389 

530 

557 

599 

2 

363 

394 

535 

0 

0 

3 

363 

389 

540 

566 

600 

4 

363 

395 

528 

0 

0 

5 

363 

0 

531 

0 

611 

6 

363 

395 

531 

564 

599 

7 

363 

390 

535 

565 

599 

8 

364 

395 

531 

557 

599 

9 

364 

390 

532 

564 

602 

10 

364 

395 

532 

564 

603 

11 

364 

395 

532 

559 

0 

12 

354 

381 

532 

563 

603 

13 

364 

395 

532 

0 

0 
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Pulse  Index 

Time  index  of  detected  peaks  for  each  multipath 

Detected 
path  #1 

Detected 
path  #2 

Detected 
path  #3 

Detected 
path  #4 

Detected 
path  #5 

14 

361 

393 

532 

562 

589 

15 

364 

396 

532 

566 

593 

16 

364 

390 

532 

0 

0 

17 

365 

396 

532 

566 

600 

18 

364 

395 

533 

0 

609 

19 

365 

391 

532 

566 

600 

20 

365 

391 

532 

559 

604 

21 

365 

396 

0 

0 

0 

22 

365 

395 

533 

0 

601 

23 

365 

396 

533 

0 

604 

24 

365 

391 

533 

564 

595 

25 

365 

396 

534 

0 

609 

26 

366 

396 

531 

565 

0 

27 

366 

396 

0 

0 

604 

28 

366 

394 

0 

0 

0 

29 

366 

392 

534 

0 

0 

Correlation 
Coefficient  ( Ef ) 

0.99999999 

0.99999997 

0.99999998 

0.99999996 

0.99999992 

Skew  ( Fj  ) 

1.00005291 

1.00004439 

1.00000411 

1.00004903 

1.00003269 

Skew  ( Fpj  ) 

-53 

-44 

-4 

-49 

-33 

Skew  error 

±17 

±30 

±21 

±44 

±61 
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Detected  Peaks  of  2  ms  LFM  chirp  from  experiment  #2  (2m) 
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Figure  3 1 .  The  peak  times  for  an  experiment  in  which  clock  skew  is  present:  The  peak 
times  ujt  derived  from  the  ensemble  of  n  =  30  measurement  signals  depicted  in 

Figure  29  are  plotted  in  this  figure.  Note  that  each  of  the  m  =  6  multipath  arrivals 

a.j ,  j  =  1, 2, . . m  has  a  corresponding  ensemble  Uj  =  j  ,  j  =  1, 2, . . . ,  m 

associated  with  it. 


Table  5.  Detailed  results  of  the  2  ms  LFM  pulse  received  during  experiment  #3,  including 
the  correlation  coefficient,  clock  skew  (in  absolute  and  ppm  values)  and  its 

associated  error. 


Analysis  of  2  ms  LFM  pulse  received  during  experiment  #3  (50  m  separation) 

Pulse  Index 

Time  index  of  detected  peaks  for  each  multipath 

Detected 
path  #1 

Detected 
path  #2 

Detected 
path  #3 

Detected 
path  #4 

Detected 
path  #5 

1 

1115 

1349 

1405 

0 

0 

2 

1115 

1357 

1406 

1450 

1523 

3 

1115 

1353 

1402 

0 

0 

4 

1116 

1358 

1401 

0 

1526 

5 

1116 

1355 

1408 

1451 

0 

6 

1116 

1351 

1401 

0 

1529 

7 

1116 

1350 

1408 

1451 

1524 

Time  index  of  detected  peaks  for  each  multipath 
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Pulse  Index 

Detected 
path  #1 

Detected 
path  #2 

Detected 
path  #3 

Detected 
path  #4 

Detected 
path  #5 

8 

1116 

0 

1408 

1451 

1524 

9 

1116 

1352 

1408 

1451 

0 

10 

1116 

0 

1413 

1443 

1524 

11 

1116 

1358 

1408 

0 

0 

12 

1116 

1351 

1410 

0 

0 

13 

1116 

1351 

1403 

1445 

1524 

14 

1116 

1352 

0 

1452 

1519 

15 

1117 

1354 

1404 

0 

0 

16 

1117 

0 

0 

1453 

0 

17 

1117 

1354 

1412 

1453 

0 

18 

1117 

1355 

1406 

0 

1524 

19 

1117 

1352 

0 

1453 

0 

20 

1117 

1352 

1405 

0 

1526 

21 

1118 

1355 

0 

1451 

1527 

22 

1118 

1352 

1409 

1449 

0 

23 

1118 

1352 

0 

1453 

1527 

24 

1118 

1352 

1412 

1453 

0 

25 

1118 

1352 

1405 

1454 

1525 

26 

1118 

0 

1412 

1448 

1529 

27 

1118 

1352 

1409 

1450 

1519 

28 

1118 

1352 

1414 

0 

1529 

29 

1118 

1357 

0 

0 

0 

Correlation 
Coefficient  ( Ej ) 

1 

0.9999987 

0.9999977 

0.9999978 

0.9999980 

Skew  ( F(  ) 

1.0000450 

0.9999948 

1.0000869 

1.0000336 

1.0000205 

Skew  ( Fpj  ) 

-45 

5 

-87 

-34 

-21 

Skew  error 

±3 

±23 

±32 

±35 

±36 
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Detected  Peaks  of  2  ms  LFM  chirp  rom  experiment  #3  (50m) 
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Figure  32.  The  peak  times  for  an  experiment  in  which  clock  skew  is  present:  The  peak 
times  ujt  derived  from  the  ensemble  of  n  =  30  measurement  signals  depicted  in 

Figure  30  are  plotted  in  this  figure.  Note  that  each  of  the  m  =  6  multipath  arrivals 

aj,j  =  \,2,...,m  has  a  corresponding  ensemble  Uj  =  jw^.  j  ,  j  =  1, 2, . . . ,  m 

associated  with  it. 


The  gradual  decrease  in  stability  of  each  successive  detected  path  is  also 
illustrated  in  Tables  4  and  5.  Unlike  the  first  path  which  is  most  likely  the  direct  path 
between  the  two  nodes,  the  following  ones  bounced  off  either  or  both  the  surface  and  the 
bottom  of  the  lake,  inducing  variability  in  the  arrivals. 

Figure  33  represents  the  received  matched-filtered  10  ms  CW  signal  recorded  in 
experiment  #3.  This  demonstrates  the  importance  of  selecting  a  probe  signal  which  can 
be  easily  distinguished  by  the  receiver.  Although  there  is  a  general  trend  in  the  shape  of 
each  response,  it  is  difficult  to  identify  a  constant  peak  over  all  of  them. 
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First  15,  10  CW  beacons  (experiment  #3) 


Figure  33.  First  fifteen  10  ms  long  CW  beacon  at  50  m  (experiment  #2). 

Finally,  Figure  34  shows  the  impact  of  selecting  a  longer  LFM  pulse.  In  this  case 
the  main  first  peak  is  still  clearly  identifiable  between  each  chirp,  but  its  general  shape  is 
larger  and  the  following  chirps  are  less  distinguishable.  The  combination  of  the  minimal 
delay  difference  between  the  arrivals  and  the  relatively  large  pulse  length  probably 
creates  too  much  superposition  to  clearly  recover  subsequent  arrivals  after  the  first  one. 
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First  15,  50ms  long  LFM  chirps  (experiment  #3) 


Time  sample  index  2000  15 


Figure  34.  First  15,  50  ms  long  LFM  chirps  at  50  m  (experiment  #3) 

Overall,  the  results  from  the  experiment  demonstrate  time  variability  in  the  tested 
channel.  The  first  arrival  is  always  present  and  strong.  However,  the  following  arrivals 
are  observed  to  be  gradually  less  stable  as  they  increasingly  interacted  with  the  channel 
boundaries.  The  algorithm  is  able  to  identify  various  multipaths  based  on  the  detected 
arrivals  and  determine  the  most  stable  one  based  on  its  correlation  coefficient.  The  time 
indices  of  these  arrivals  clearly  illustrate  the  presence  of  a  skew  between  the  two  clocks 
as  the  time  index  of  the  main  arrival  slowly  shifts  to  the  right,  as  shown  in  Tables  4  and  5 
and  the  associated  Figures  3 1  and  32.  The  variance  of  the  results  in  Table  3  could  be  due 
the  short  period  of  time  over  which  the  chirps  are  transmitted. 
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VII.  CONCLUSION  AND  RECOMMENDATIONS 


In  this  research  we  addressed  the  problem  of  clock  synchronization  in  an 
underwater  acoustic  network.  This  is  a  very  actively  researched  area  for  every  type  of 
communication  channel.  This  is  particularly  challenging  in  the  underwater  acoustic 
environment  due  to  its  high  latency  characteristics,  and  it  requires  the  development  of  a 
specific  clock-synchronization  protocol.  Most  of  the  approaches  found  in  the  literature 
assumed  constant  propagation  delays  between  two  nodes  and  ignore  time  variability  of 
the  channel.  An  overview  of  how  and  why  time  variance  exists  in  the  underwater 
acoustic  environment  was  given  and  a  method  to  evaluate  and  mitigate  its  effects  as  part 
of  a  clock-synchronization  protocol  was  proposed. 

The  proposed  method  is  divided  in  two  stages.  The  first  stage  evaluates  the 
frequency  drift  between  the  clocks  of  two  adjacent  nodes  (called  the  “skew”),  and  the 
second  stage  estimates  the  offset  between  the  same  two  nodes. 

The  first  stage  was  simulated  using  a  computer  model  and  tested  in  a  small, 
shallow,  fresh  water  lake.  Both  simulated  and  experimental  results  show  the 
effectiveness  of  the  proposed  method.  However,  further  testing  needs  to  be  done  on  the 
basis  of  the  following  recommendations  to  address  the  variability  in  the  estimated  skew 
calculated  from  the  experimental  results. 

The  second  stage  could  only  be  simulated  since  an  experiment  would  require 
modifying  the  firmware  of  the  modem,  which  was  not  readily  possible  with  the  Teledyne 
Benthos  acoustic  modems.  The  results  from  the  simulations  are  in  accordance  with  the 
predictions  from  the  theoretical  model  of  the  proposed  method  for  offset  correction. 

Additional  experiments  are  required  to  fully  develop  the  method  and  select  the 
optimal  parameters.  Once  completely  tested,  the  addition  of  this  capability  to  an 
underwater  acoustic  network  will  keep  the  clocks  of  the  nodes  synchronized.  This 
synchronization  should  improve  the  utility  of  the  network  when  deployed  or  used  for 
long  periods  of  time. 
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The  following  recommendations  do  not  necessarily  result  from  the  simulation  and 
experiments.  Some  of  them  are  simply  aspects  that  were  omitted  due  to  time  or  resource 
constraints. 

A.  REPRODUCE  EXPERIMENT  FROM  LITERATURE 

Most  of  the  literature  pertaining  to  underwater  clock  synchronization  does  not 
describe  physical  experiments  in  an  actual  underwater  environment.  The  simulations 
reported  in  [1]  and  [7]  should  be  reproduced  using  the  protocol  explained  in  this  thesis. 
This  would  allow  for  some  comparison  among  the  various  protocols  and  confirm  that 
accounting  for  time  variability  improves  the  synchronization  of  the  clocks. 

Another  approach  to  compare  the  different  protocols  is  to  reproduce  our 
experiment  using  the  protocols  from  [1]  and  [7]. 

B.  PSEUDORANDOM  BINARY  SEQUENCE  (PRBS) 

Another  type  of  pulse  used  in  channel  sounding  is  the  pseudorandom  binary 
sequence  (PRBS)  [10].  Just  as  the  LFM  chirp,  the  PRBS  is  characterized  by  its  narrow 
ambiguity  function,  making  it  an  excellent  choice  for  measuring  the  impulse  response  of 
a  channel.  Its  structure  allows  it  to  be  more  tolerant  to  superposition.  The  generation  of  a 
PRBS  signal  is  not  complicated  and  reproducing  the  experiment  from  Chapter  VI  would 
require  minimal  effort. 

C.  LONGER  RPI 

As  discussed  in  Chapter  VI,  the  variability  of  the  results  shown  in  Table  3  might 
originate  from  a  short  RPI.  The  experiment  should  be  reproduced  using  a  constant  2  ms 
LFM  chirp  while  increasing  the  RPI  duration.  This  would  confirm  our  hypothesis  that  the 
results  would  improve  with  a  longer  RPI. 

D.  OFFSET  CORRECTION  EXPERIMENT 

The  simulation  proves  the  theory  of  the  offset  correction  algorithm.  However,  an 
experiment  in  the  underwater  environment  can  possibly  highlight  any  neglected  factors  as 
was  shown  when  testing  our  skew-correction  algorithms.  As  stated  earlier,  technical  and 
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time  limitations  prevented  us  from  testing  our  offset  correction  algorithm  in  Del  Monte 
Lake.  An  experiment  should  be  designed  and  conducted  to  test  the  developed  algorithm 
for  offset  correction  described  in  Chapter  IV  and  simulated  in  Chapter  V. 

E.  MORE  REALISTIC  ENVIRONMENT 

Del  Monte  Lake  was  selected  mostly  for  its  accessibility.  As  stated  in  Chapter  III, 
the  effects  of  time  variability  are  location  dependent.  Since  a  shallow  fresh  water  lake  is 
not  a  probable  environment  where  an  underwater  acoustic  network  would  be  deployed, 
the  experiment  should  be  reproduced  in  an  environment  more  typical  of  where  a  time- 
synchronization  protocol  would  be  required  and  beneficial. 

One  of  the  advantages  of  Del  Monte  Lake  was  its  low  noise  level.  This  will  not 
necessarily  be  true  in  an  ocean  environment.  Future  work  should  involve  folding 
measurement  noise  into  the  model  for  skew  and  offset.  Also,  noise  reduction  algorithms 
should  be  explored. 

F.  DETECTED  ARRIVALS  SORTING  ALGORITHM 

The  current  algorithm  is  vulnerable  to  noise  and  superimposed  arrivals.  Future 
work  should  explore  new  algorithms  to  cope  with  these  issues. 
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